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Objectives: 

The  objectives  of  this  project  have  not  changed  since  the  last  status  report. 

Status  of  Effort: 

This  Phase  II  STTR  project  was  completed  on  schedule  on  September  14,  2001. 

Accomplishments  /  New  Findings: 

This  report  covers  the  period  between  June  15,  1999  and  September  14,  2001.  Detailed 
technical  information  on  the  research  described  is  not  presented  here,  but  can  be  found  in 
attached  reports  and  journal  publications  and  in  portions  of  the  recent  text  edited  by  R.  J.  Barker 
and  E.  Schamiloglu  [1].  Many  of  the  algorithms  discussed  here  have  been  implemented  in  the 
XOOPIC  code  [2],  which  is  freely  available  for  noncommercial  use  -  thus,  providing  benefit  to 
the  entire  PIC  community. 


The  Electron  Self-Force 

In  this  section,  work  on  the  electron  self-force  issue  is  summarized.  The  attached  report,  “A 
Radiation  Damping  Algorithm  for  Particle  Simulation”,  describes  this  work  in  more  detail.  This 
report  will  be  submitted  for  publication  by  the  authors. 

The  electron  self-force  cannot  be  modeled  by  a  conventional  PIC  algorithm,  due  to  the  finite 
resolution  of  temporal  and  spatial  scales.  A  sophisticated  analytical  treatment  of  self-force 
effects  on  charged  particles  has  been  developed  previously  [3].  However,  a  recent  report  [4]  has 
applied  this  work  to  the  problem  of  an  electron  beam  in  an  HPM  tube,  showing  that  self-force 
effects  are  negligible  in  microwave  HPM  devices. 

The  full  equations  from  Ref.  [3]  describing  the  self-force  are  quite  complicated  and  would  be 
computationally  expensive  to  include  in  a  PIC  algorithm.  More  seriously,  these  equations 
include  the  second  time  derivative  of  the  particle  momentum  --  a  very  noisy  quantity  in  a  time- 
explicit  electromagnetic  PIC  algorithm.  Thus,  it  is  not  practical  to  directly  implement  the  results 
of  Ref.  [3]  in  a  PIC  code. 

In  the  physical  limit  where  self-force  effects  are  important  (i.e.,  for  ultra-relativistic 
electrons),  one  can  approximate  the  self-force  as  the  classical  recoil  of  the  particle  upon  emitting 
synchrotron  radiation.  The  equations  describing  this  effect  [5]  are  relatively  simple  and  require 
only  the  first  time  derivative  of  the  particle  momentum.  These  equations  have  been  successfully 
implemented  and  tested  in  XOOPIC.  It  is  unlikely  that  this  work  will  improve  PIC  simulations 
of  HPM  tubes  in  the  near  term;  however,  it  is  relevant  to  beam  devices  in  which  ultra-relativistic 
electrons  interact  with  strong  magnetic  fields. 
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We  also  explored  the  issue  of  synchrotron  radiation  effects  (i.e.  bremsstrahlung)  during 
collisions  between  beam  electrons  and  a  tenuous  background  neutral  gas.  An  unpublished  report 
[6]  directly  addresses  this  issue  for  many  monatomic  gasses.  This  physics  could  be  accurately 
modeled  in  XOOPIC  through  a  Monte  Carlo  collision  mechanism,  which  would  be  more 
appropriate  than  any  direct  modification  to  the  electromagnetic  PIC  algorithm.  However,  Ref. 
[6]  shows  that  energy  losses  due  to  bremsstrahlung  are  negligible  until  the  electron  energy 
exceeds  100  MeV  and  the  neutral  gas  density  becomes  much  larger  than  is  found  in  an  evacuated 
tube.  Thus,  we  do  not  plan  to  further  pursue  this  issue. 

Secondary  Emission  of  Electrons 

In  this  section,  the  work  on  the  secondary  emission  model  is  summarized.  The  complete 
details  of  the  model  can  be  found  in  Sec.  11.4.2,  beginning  on  p.  397  of  Ref.  [1], 

Secondary  emission  can  play  an  important  role  in  many  devices  where  electrons,  ions,  or 
neutral  particles  impact  surfaces.  The  impact  can  result  in  reflection  of  the  incident  particle, 
multiple  scattering  events  within  the  surface  lattice  followed  by  ejection,  or  transfer  of  the 
incident  energy  to  electrons  in  the  lattice  leading  to  ejection  of  one  or  more  electrons. 

The  secondary  model,  implemented  in  the  XOOPIC  code,  includes  a  simple  model  and  a 
more  complete  model.  In  the  simple  model,  an  energy  and  angular-independent  secondary 
emission  coefficient  is  specified.  The  emitted  electrons  are  modeled  as  true  secondaries,  ejected 
with  a  uniform  angular  distribution  at  a  specified  temperature.  The  complete  model  allows 
specification  of  the  key  parameters  to  model  energy  and  angular  dependence,  relevant  to  most 
conductive  materials  and  many  dielectric  materials  as  well.  The  complete  model  also  considers 
surface  roughness  in  computing  the  secondary  emission  coefficient.  In  the  complete  model,  the 
fractions  of  reflected  primaries,  scattered  primaries,  and  true  secondaries  can  be  specified,  along 
with  the  energy  and  angular  distribution  functions  for  the  true  secondaries.  The  current  model 
defaults  to  isotropic  angular  distribution,  which  is  appropriate  for  most  metals.  Fractional  yields 
are  emitted  statistically  using  pseudorandom  numbers. 

Secondary  emission  plays  an  important  role  in  many  devices  relevant  to  microwave  sources, 
including  two-surface  RF  multipactor  phenomena,  single  surface  DC  multipactor  phenomena 
(such  as  window  breakdown),  as  well  as  beam  intercept  current  in  collectors  or  circuit  structures 
in  microwave  sources  and  accelerators,  and  ion/neutral-induced  secondary  emission  in  DC  and 
RF  discharges. 


Loading  and  Injection  Models 

This  section  summarizes  the  work  on  loading  and  injection.  A  careful  analysis  of  the 
methods  for  loading  and  injecting  particles  in  PIC  codes  has  indicated  deficiencies  in  the 
accuracy  of  common  methods.  This  work  is  described  in  detail  in  Ref.  [7],  where  the  analysis  of 
existing  methods  is  performed,  and  correction  terms  are  proposed  for  second  order  accuracy  for 
loading  and  injection  in  the  general  case. 

A  second-order  accurate  method  for  loading  general  phase  space  distributions  of  particles  has 
been  described.  The  method  for  time  centering  described  in  Birdsall  and  Landon  [8]  is  shown  to 
be  first  order  accurate  when  the  time-derivative  of  the  acceleration  is  non-zero,  the  magnetic  field 
is  non-zero,  or  the  initial  velocity  is  non-zero.  A  correction  is  derived  to  make  ensure  second 
order  accurate  time  centering  of  loaded  particles  for  the  general  case. 

The  injection  of  particle  flux  from  a  boundary  in  the  system  was  also  considered.  The 
injection  may  be  due  to  a  thermionic  or  field  emitter,  or  it  may  be  flux  crossing  a  mathematical 
boundary  line  with  a  specified  distribution  outside  the  system.  Analysis  indicates  that  existing 
injection  schemes  result  in  a  first  order  error  term.  In  this  work,  second  order  accuracy  was 
achieved  for  a  number  of  special  cases  by  adding  additional  terms  to  cancel  the  first  order  error. 
It  was  found  that  cancellation  of  error  terms  leads  to  an  increasingly  complicated  algorithm  as  the 
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complexity  of  the  terms  and  their  derivatives  in  time  and  space  increases.  Analysis  of  the 
computational  efficiency  of  the  various  correction  schemes  was  also  performed. 

The  loading  and  injection  schemes  are  relevant  to  applications  in  which  particles  are  loaded 
initially  (such  as  discharges  or  plasma-filled  microwave  tubes),  or  injected  from  the  boundaries 
of  the  system  (such  as  the  cathode  in  a  microwave  tube,  or  secondary  emission  from  a  collector). 
A  number  of  corrective  schemes  are  presented  which  minimize  the  increased  complexity  of  the 
equations  of  motion  for  achieving  second  order  accuracy. 

Accuracy  Analysis  of  the  PIC  Method 

In  this  section,  the  analysis  of  the  accuracy  of  the  PIC  model  is  summarized.  More  details  are 
provided  in  three  presentations  [9,10,1 1],  all  of  which  are  attached  to  this  report. 

The  discretization  errors  for  each  component  stage  of  the  PIC  method  are  derived  for  both 
electrostatic  and  electromagnetic  models.  Some  of  the  error  terms  have  been  described 
previously;  here  the  analysis  is  extended  to  include  errors  in  the  interpolation  algorithm,  which 
connect  the  continuum  particles  to  the  discrete  fields.  The  terms  coupling  errors  between  PIC 
steps  are  also  described.  The  error  analysis  includes  both  spatial  and  temporal  errors. 

For  the  electrostatic  PIC  model,  the  steps  analyzed  include  the  interpolation  of  position  to  a 
mesh  to  compute  the  charge  density,  the  discretization  error  of  the  finite  difference  solution  of 
Poisson’s  equation,  errors  in  computing  the  gridded  electric  field  from  potential,  errors  in 
interpolating  the  electric  field  from  the  mesh  to  the  continuum  particle  locations,  and  the  error 
terms  in  integrating  Lorentz’s  equation.  Finally  the  individual  errors  are  incorporated  into  a 
coupled  error  term.  It  was  particularly  interesting  to  note  that  the  error  terms  in  the  linear  charge 
weighting  precisely  cancel  those  in  the  Poisson  equation.  The  total  error  for  the  linear  weighting 
case  was  shown  to  contain  error  terms  no  worse  than  second  order  in  both  space  and  time,  with 
dependence  upon  higher  derivatives  of  the  potential. 

For  the  electromagnetic  PIC  model,  the  steps  analyzed  include  the  interpolation  from  the 
continuum  positions  and  velocities  of  particles  to  the  discrete  locations  of  the  Yee  mesh  to 
compute  the  current  density,  the  errors  in  the  finite  difference  discretization  of  Ampere’s  law, 
and  Faraday’s  law,  the  errors  in  the  interpolation  from  the  Yee  mesh  location  to  common  nodal 
locations,  the  errors  from  interpolating  the  nodal  fields  to  the  continuum  particle  positions,  and 
the  integration  of  the  relativistic  Loren tz  equation. 

Charge  Weighting  Correction 

This  section  provides  a  summary  of  the  work  on  charge  weighting  corrections  for  non- 
cartesian  geometry.  Complete  details  are  provided  in  Ref.  [12],  which  is  attached  to  this  report. 

This  novel  formulation  of  the  charge  weighting  algorithm  for  the  particle-in-cell  method  was 
developed  to  eliminate  a  systematic  error  in  the  charge  density  in  curvilinear  systems  which 
occurs  for  the  classic  PIC  weighting  method  [8]. 

In  the  standard  PIC  scheme,  particles  exist  in  a  continuum  phase  space,  while  the  fields  are 
defined  on  a  discrete  mesh.  The  particles  and  fields  interact  through  various  interpolation 
schemes.  The  source  terms  for  the  field  equations  are  the  current  density,  J,  and  the  charge 
density,  p.  These  quantities  are  computed  weighting  the  current  or  charge  for  each  particle  to 
mesh  nodes,  and  then  dividing  by  the  appropriate  geometric  terms  to  obtain  the  current  or  charge 
density.  The  geometric  terms  were  viewed  as  a  property  of  the  mesh,  independent  of  the 
interpolation  scheme  used  for  the  particles.  These  schemes  lead  to  systematic  errors  of  33%  and 
100%  on  axis  in  cylindrical  and  spherical  coordinates,  respectively,  which  remain  even  as  the 
mesh  size  approaches  the  infinitesimal  limit.  Errors  also  occur  at  the  outer  radial  edge,  and  at 
interior  boundaries  as  well.  The  errors  occur  at  every  mesh  point  for  non-uniform  meshes.  Each 
of  these  errors  is  shown  in  Figs.  1  and  2  in  the  attached  preprint  [12],  Case-dependent  correction 
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factors  have  been  described  previously,  but  each  of  these  is  valid  only  for  a  specific  mesh  with  a 
specific  interpolation  scheme. 

In  this  work,  we  derive  an  improved  weighting  algorithm  which  requires  no  correction  for  the 
general  case,  valid  for  any  mesh  with  any  interpolation  scheme.  Rather  than  computing  the 
volume  or  surface  area  elements  from  purely  geometric  considerations,  this  scheme  employs  an 
interpolation  for  computing  the  volume  elements.  The  interpolation  is  symmetric  with  that  used 
for  weighting  the  charge/current  to  the  mesh,  and  removes  all  systematic  errors.  The  method  has 
the  desirable  properties  that  total  volume  is  conserved,  as  is  total  charge/current,  and  it  yields  the 
exact  answer  for  the  uniform  charge  distribution.  Arbitrary  distributions  are  correct  to  the  mesh 
resolution,  with  error  approaching  zero  as  the  mesh  size  approaches  zero. 

This  method  is  relevant  to  any  plasma  device  model  in  curvilinear  coordinates,  such  as  a 
cylindrical  model  of  an  O-type  microwave  tube.  The  new  scheme  will  result  in  reduced  noise  and 
errors  for  particles  near  the  axis  and  near  any  boundary  as  well,  which  means  a  less  noisy  RF 
output  signal  in  the  above  case. 

Digital  Filter  Noise  Reduction 

This  section  provides  a  summary  of  the  work  on  digital  filter  noise  reduction.  Complete 
details  are  provided  in  Ref.  [13],  which  is  attached  to  this  report.. 

Digital  filters  are  used  in  many  digital  applications  to  reduce  high  frequency  or  short 
wavelength  noise.  The  technique  has  been  described  for  use  in  particle  simulation  for  Cartesian 
coordinates  [Birdsall85].  However,  the  method  cannot  be  applied  directly  to  particle  simulation 
in  cylindrical  coordinates,  a  domain  of  great  interest  in  microwave-beam  devices,  beam  optics, 
plasma-filled  devices,  and  many  others.  In  these  cylindrical  systems,  the  noise  on  axis  can  be 
considerable  when  the  number  of  simulation  particles  per  cell  is  small;  this  is  exacerbated  by  the 
cylindrical  coordinate  system  where  the  number  of  particles  per  cell  is  proportional  to  the  radius 
for  a  uniform  density  or  current  profile.  In  addition  to  the  immediate  effects  of  noise,  the  short 
wavelength  noise  associated  with  low  numbers  of  computer  particles  per  cell  can  lead  to 
numerical  heating. 

We  have  extended  the  standard  1-2-1  digital  filter  to  the  cylindrical  coordinate  system,  and 
the  model  has  been  implemented  and  tested  in  one  dimension;  it  will  be  extended  to  two 
dimensions  as  part  of  an  ongoing  study.  The  filter  operates  on  charge  density  in  an  electrostatic 
model,  and  on  current  density  in  an  electromagnetic  model. 

The  filter  was  designed  to  satisfy  a  number  of  desirable  criteria.  It  conserves  net  charge  and 
current  exactly.  It  preserves  a  uniform  distribution;  i.e.  charge/current  does  not  migrate 
systematically  towards  the  axis  or  outer  radius.  Furthermore,  analysis  indicates  the  filter  has 
strength  decreasing  rapidly  with  radius,  which  can  provide  the  benefits  of  filtering  close  to  the 
axis  with  minimum  impact  of  the  sheath  or  beam  profile  at  the  outer  radius.  Due  to  the  weakness 
of  the  filter,  it  was  empirically  observed  that  hundreds  of  passes  were  required  for  reduction  of 
short  wavelength  noise  that  five  to  10  passes  of  a  planar  1-2-1  filter  would  provide.  A  pre¬ 
computed  multi-pass  filter  matrix  was  devised  which  would  apply  an  arbitrary  number  of  passes 
with  a  single  matrix  multiply.  The  multi-pass  filter  was  demonstrated  to  be  far  effective  at 
reducing  noise  than  the  addition  of  particles. 

The  multi-pass  digital  filter  described  in  the  attached  manuscript  provides  a  number  of 
benefits  relevant  to  microwave-beam  devices.  It  can  be  used  to  reduce  particle-generated  noise  to 
acceptable  levels,  as  well  as  reduction  of  noise  and  numerical-heating  effects  in  plasmas  in 
devices  such  as  plasma-filled  or  focused  microwave  tubes,  or  in  plasma  discharges  in  cylindrical 
configurations. 
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Fluctuation  Reduction  Algorithms 

This  section  provides  a  summary  of  the  work  on  fluctuation  reduction  algorithms.  Complete 
details  are  provided  in  the  final  report  from  U.  Michigan  [14]  and  the  presentation  by  Prof.  Boyd 
[15],  both  of  which  are  attached  to  this  report.. 

The  objective  of  this  work  was  to  demonstrate  the  application  of  a  direct  simulation  Monte 
Carlo  (DSMC)  [16]  algorithm,  called  information  preservation  (IP),  [17,18]  to  noise  reduction  in 
a  1-D  particle-in-cell  (PIC)  code.  This  technique  was  successfully  implemented  in  the  1-D 
electrostatic  PIC  code  XES1,  and  effective  noise  reduction  has  been  shown. 

A  new  Graphical  User  Interface  for  XOOPIC 

The  initial  version  of  the  new  cross-platform  graphical  user  interface  (GUI)  for  XOOPIC  is 
complete.  This  professional  GUI  runs  equally  well  under  Microsoft  Windows,  Linux  or  any 
version  of  Unix.  Work  on  the  GUI  is  continuing  under  a  project  funded  by  DoE,  and  we  plan  to 
support  the  Macintosh  OSX  in  the  near  future. 

The  GUI  provides  users  with  full  control  of  the  simulation  and  displays  real-time  2-D  and 
3-D  plots  of  relevant  physical  quantities.  This  interactive  approach  to  PIC  simulation  provides 
insights  that  are  difficult  to  obtain  from  the  post-processing  of  data. 

The  commercial  version  of  XOOPIC,  to  be  called  OOPIC  Pro,  is  now  ready  for  beta  testing. 
A  commercial  licensing  agreement  has  been  negotiated  with  the  University  of  California 
Berkeley,  enabling  Tech-X  to  satisfy  the  commercialization  and  technology  transfer  aspects  of 
the  STTR  program. 

An  initial  user  manual  was  developed  for  OOPIC  Pro  under  this  contract.  Further  work  on 
the  user  manual  is  continuing  under  a  project  funded  by  DoE.  An  excerpt  from  the  user  manual 
is  provided  here,  in  order  to  show  what  has  been  accomplished. 

Excerpt  from  the  OOPIC  Pro  user  manual: 

Windows  users  begin  using  the  code  by  double-clicking  on  the  oopicpro.exe  executable  file.  The 
following  file  dialog  will  appear: 
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Browse  until  you  have  found  the  directory  containing  the  input  file  you  wish  to  open.  The  defualt 
directory  is  oopicpro/input.  Double  click  on  the  desired  input  file.  Alternatively,  click  once  on 
the  file  you  want  to  open  and  click  Open,  or  type  the  name  of  the  file  and  hit  <RETURN>.  Upon 
choosing  an  input  file,  the  main  OOPIC  Pro  control  window  will  appear: 
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At  the  top  of  the  control  window,  there  is  a  menu  bar.  Below  the  bar  are  four  buttons  which 
control  the  progress  of  the  simulation.  In  the  interior  of  the  window  are  the  two  program  logs. 
The  current  simulation  time  is  given  on  the  left  side  of  the  lower  window  border,  and  the  name  of 
the  current  input  file  on  the  right  side  of  the  border. 

File  menu 

Clicking  on  the  File  menu  gives  the  user  access  to  Loading  and  Saving  files.  Although  the  user 
must  specify  an  input  file  upon  starting  the  program,  one  can  Open  another  input  file  in  order  to 
begin  a  new  simulation. 

The  File  menu  also  contains  an  Exit  heading.  Clicking  on  Exit  brings  up  a  query  window  asking 
the  user  if  he  wants  to  end  the  program.  One  can  also  end  the  program  immediately  by  clicking 
on  the  ‘X’  button  in  the  upper  right  hand  corner  of  the  control  window. 

Edit,  menu 

The  Edit  menu  allows  the  user  to  Cut,  Copy  and  Paste.  The  user  can  also  view  the  input  file  in  a 
read-only  browser  window  by  clicking  on  Input  text  file  under  the  Edit  window. 

Under  the  View  menu,  choosing  Style  brings  up  a  sub-menu  with  choices  about  the  appearance 
of  the  windows  in  OOPIC  Pro. 

Also  under  the  Vi  ew->Logs  submenu,  one  can  choose  whether  to  view  the  program  logs. 
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After  choosing  the  last  heading  under  the  View  menu,  Diagnostics,  the  list  of  available 
diagnostic  plots  is  displayed.  Clicking  on  any  of  these  will  bring  up  the  corresponding  diagnostic 
plot. 

Parameter  menu. 

The  first  sub-menu  under  the  parameter  menu  is  denoted  Iterations.  Clicking  on  Maximum 
Number  in  that  sub-menu  brings  up  an  value  entry  window,  in  which  the  user  is  able  to  set  the 
maximum  number  of  time  steps  in  the  simulation.  The  default  value  is  ‘0‘,  which  specifies  that 
there  is  no  limit  on  the  number  of  iterations.  Change  the  value  by  typing  another  integer,  or  by 
clicking  on  the  up  or  down  arrows  on  the  right  edge  of  the  value  blank.  When  the  entry  blank 
contains  the  desired  number,  click  the  Ok  button  or  hit  <RETURN>. 

Choosing  Periodic  dump  brings  up  a  similar  value  entry  window  in  which  the  user  can  choose 
how  often  to  periodically  save  the  progress  of  the  simulation.  The  default  is  ‘O’,  meaning  no 
dump  files  are  saved.  More  information  about  dump  files  later  in  this  section  of  the  manual. 

Clicking  on  Number  per  update  causes  another  value  entry  window  to  appear,  in  which  the  user 
can  change  the  how  many  time  steps  are  between  graphics  updates. 

The  second  sub-menu  in  the  Parameter  menu  is  titled  Movies.  These  parameters,  as  well  as  the 
entire  process  for  creating  a  movie  of  the  plots  in  a  simulation,  will  be  explained  later  in  this 
section. 

The  next  and  last  menu,  the  Help  menu,  allows  the  user  to  view  various  documentation  for  the 
program.  One  can  access  this  manual  by  choosing  Manual.  Clicking  on  Software  Documentation 
brings  up  code-specific  documentation  produced  by  Doxygen.  Clicking  the  Text  Documentation 
header  causes  the  doc  directory  window  to  appear.  It  contains  text  files  which  were  the 
documentation  for  OOPIC  and  XOOPIC,  which  can  be  viewed  using  any  text  viewer.  All 
currently  applicable  information  from  the  text  files  have  been  incorporated  into  this  document. 
Copyright  information  can  be  viewed  by  clicking  on  the  About  heading. 

VCR  Control  Buttons 

The  four  large  VCR  buttons  underneath  the  menu  bar  control  the  execution  of  the  simulation. 
The  [Run]  button  causes  the  simulation  to  run.  The  [Step]  button  advances  the  simulation  one 
frame  (note  that  this  may  not  be  the  same  as  advancing  the  simulation  one  time  step.  The  [Stop] 
button  stops  the  simulation.  The  [Restart]  button  resets  the  simulation  back  to  the  initial 
conditions. 

The  diagnostic  plots: 

Upon  opening  OOPIC  Pro,  there  will  be  position  plots  for  each  of  the  species  present  in  the 
simulation.  The  user  can  bring  up  plots  of  other  diagnostic  parameters.  First,  open  the  View 
menu  (click  on  Vi  ew  on  the  menu  bar),  hold  the  mouse  arrow  briefly  over  the  word  Diagnostic 
which  appears  with  the  rest  of  the  menu  (see  the  screen  shot  below).  Then  point  the  mouse  arrow 
to  the  desired  diagnostic  and  click  again.  A  new  window  containing  the  plot  will  appear.  Note: 
this  can  be  done  while  the  simulation  is  running.  If  the  same  diagnostic  is  selected  a  second 
time,  the  plot  is  removed. 
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There  are  four  types  of  dignostic  plots  -  particle  plots,  3-D  surface  plots,  line  plots  and  vector 
plots. 
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Particle  plots. 

Particle  plots  include  the  configuration  plots,  which  show  the  2-D  simulations  region,  with 
physical  boundaries  shown  in  various  colors  and  the  macroparticles  of  a  given  species  shown  as 
colored  points.  There  is  one  configuration  plot  for  each  particle  species  -  an  example  is  provided 
by  the  screen  shot  below.  There  are  also  phase  space  particle  plots,  which  show  the 
macroparticle  locations  with  a  spatial  coordinate  along  the  horizontal  axis  and  a  momentum 
coordinate  along  the  vertical  axis.  One  plot  for  each  combination  of  position  and  momentum  is 
provided.  In  the  phase  space  plots,  the  macroparticles  for  all  species  are  shown  simultaneously, 
differentiated  by  color.  These  colors  are  used  consistently  in  all  of  the  particle  plots. 


Like  all  of  the  2-D  plots  in  the  OOPIC  Pro  GUI,  the  user  can  interactive  change  the  axes  limits, 
axes  labels,  plot  title,  etc.  One  can  also  zoom  in  to  a  smaller  area  of  the  plot,  by  choosing  Zoom 
i  n  from  the  Vi  ew  menu  and  then  selecting  the  desired  region  with  the  mouse.  After  multiple 
zoom  operations,  the  user  can  scroll  back  through  previous  views  by  clicking  the  left  mouse 
button.  At  any  point,  the  user  can  choose  to  apply  the  current  settings  and  continue  with  the 
simulation. 

Line  plots. 

OOPIC  Pro  also  provides  a  number  of  line  plots,  which  show  the  time  history  of  a  specified 
quantity  -  an  example  is  provided  by  the  screen  shot  below,  which  shows  the  number  of 
macroparticles  in  the  simulation. 
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3-D  surface  plots.  « 

OOPIC  Pro  provides  many  3-D  surface  plots  of  quantities  defined  on  the  2-Dgn-® 
is  provided  by  the  screen  shot  below,  which  shows  the  component  f  h  1  neutral 
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Vector  plots. 

OOPIC  Pro  also  provides  2-D  vector  plots  of  quantities  defined  on  the  2-D  grid  -  an  example  is 
provided  by  the  screen  shot  below,  which  shows  the  radial  component  of  the  electric  field.  In 
addition  to  all  of  the  field  components,  one  can  plot  particle  densities,  the  density  of  any  neutral 
gasses  in  the  simulation,  or  FFT  power  spectral  densities. 
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Saving  and  printing  plots 

Saving  and  printing  are  the  same  for  all  types  of  plots.  To  save  a  plot,  select  Save  as  from  the 
File  menu.  In  the  window  which  appears,  type  a  file  name  and  choose  a  format  in  which  to  save 
the  image.  To  print  a  plot,  choose  P  ri  nt  from  the  File  menu. 

Log  windows 

Both  the  Standard  output  log  and  the  Standard  error  log  are  present  at  startup  by  default.  The 
standard  output  log  is  in  the  upper  part  of  the  interior  of  the  main  control  window.  It  displays 
comments  on  the  status  of  the  program,  and  gives  messages  as  each  component  of  an  input  file  is 
interpreted  and  added  to  the  simulation.  The  standard  error  file  is  in  the  lower  part  of  the  interior 
of  the  main  control  window. 

Restart  files 

One  can  save  the  current  state  of  a  simulation  in  a  cross-platform  binary  restart  file.  A  restart  file 
is  in  a  block  format  where  the  specifications  of  the  simulation  and  the  current  values  of  particle 
positions,  fields,  boundaries,  and  diagnostic  values  such  as  time  histories,  are  written.  To  save  a 
restart  file,  click  on  the  File  menu,  choose  the  Save  header,  then  click  on  Dump  file.  The  default 
name  for  the  dump  file  is  the  input  file  name  with  a  .dmp  extension. 

One  can  also  set  OOPIC  Pro  to  periodically  save  a  restart  file  during  a  simulation.  Under  the 
Parameter  menu,  choose  Iterations  and  click  on  Periodic  Dump.  Then  change  the  number  to  be 
how  many  time  steps  are  desired  between  file  saves.  The  user  can  specify  the  name  of  the  file. 
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To  restore  a  simulation  to  the  state  at  which  it  was  when  the  restart  file  was  saved,  choose  Open 
under  the  File  menu  and  click  on  Dump  file.  Then  open  the  directory  containing  the  desired 
dump  file  and  double  click  on  the  file.  Note  that  although  the  simulation  will  be  at  the  same  state 
as  before,  GUI  will  not  return  to  the  previous  state,  since  no  information  about  the  GUI  is 
recorded  in  the  dump  file.  Note:  Currently,  after  loading  a  dump  file,  if  one  wants  to  completely 
restart  with  the  initial  conditions  of  the  simulation,  or  to  load  another  input  file,  OOPIC  Pro  must 
be  exited  and  re-run. 

The  movie  feature: 

OOPIC  Pro  has  the  capability  of  saving  multiple  frames  of  one  or  more  diagnostic  plots  during  a 
simulation  for  later  viewing  as  a  movie.  The  following  is  a  description  of  the  process  of  making 
a  movie. 

First,  in  the  Parameters  menu,  choose  Movies,  and  then  click  on  “Save  Images  as”.  By  default, 
OOPIC  Pro  will  place  the  image  files  in  a  directory  named  the  same  as  the  input  file.  If  one  is 
making  a  movie  of  a  simulation  for  the  first  time,  a  box  will  pop  up,  informing  the  user  that  a 
new  directory  called  cinput  file>_movie  is  being  created,  wherecinput  file>  is  the  name  of  the 
input  file  without  the  .inp  extension.  Hit  OK,  and  a  dialog  window  will  appear.  The  default 
image  file  stub  appears  in  the  File  Name  input  line.  It  is  the  input  file  name  followed  by  an 
underscore.  Change  the  file  stub  if  desired.  Then  choose  the  format  of  the  movie  image  files. 
The  actual  image  file  names  will  consist  of  the  image  file  stub,  followed  by  the  name  of  the  plot 
being  saved,  followed  by  the  frame  number,  with  an  image  file  extension.  For  example,  a  typical 
file  name  would  be  “cavity_z-r_phase_space_for_electrons.00001.jpeg”. 

If  the  image  file  stub  is  not  set  by  the  user,  the  default  is  the  input  file  name.  The  default 
directory  is  the  oopicpro/input  directory  (not  oopicpro/input/cinput  file  name>_movie,  and  the 
default  image  format  is  ‘.ppm’. 

After  choosing  an  image  file  stub,  the  user  must  choose  the  frame  rate  of  the  movie.  Again  under 
the  Parameters  menu  choose  the  Movies  sub-menu,  and  click  on  Iterations  per 
Output.  The  number  in  the  entry  window  which  appears  specifies  the  number  of  simulation 
plot  frames  per  movie  frame.  In  a  similar  fashion,  one  can  also  specify  the  maximum  number  of 
frames  in  the  movie. 

If  the  simulation  has  multiple  graphic  plots  open,  a  set  of  image  files  will  be  created  from  each 
plot  window.  Rather  than  load  all  of  these  files  into  the  movie  displaying  program  sequentially, 
the  user  can  specify  a  Config  file,  which  keeps  track  of  the  various  plot  movies. 

Personnel  Supported: 

This  project  supported  personnel  from  three  institutions,  including  Tech-X  Corporation,  the 
Plasma  Theory  and  Simulation  Group  (PTSG)  at  the  University  of  California  Berkeley,  and  a 
group  at  the  University  of  Michigan. 

Supported  personnel  at  Tech-X  included  the  principal  investigator  Dr.  David  Bruhwiler,  as 
well  as  Kelly  Luetkemeyer  and  Dr.  John  Cary. 

Supported  personnel  at  UC  Berkeley  included  the  principal  investigator  Prof.  C.  K.  Birdsall, 
as  well  as  Dr.  John  Verboncoeur  and  students  Xingbo  Yu  and  Peter  Mardahl.  Dr.  Keith 
Cartwright  completed  his  PhD  thesis  work  with  partial  support  from  this  project,  and  is  now  at 
the  Air  Force  Research  Laboratory  (AFRL)  at  Kirtland  AFB  in  New  Mexico. 
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Supported  personnel  at  the  University  of  Michigan  included  the  principal  investigator  Prof. 
Iain  Boyd,  as  well  as  Dr.  Jing  Fan. 

Publications: 

See  the  “attachments”  subsection  below.  Also,  much  of  the  work  funded  under  this  project 
has  appeared  in  the  recent  text  “High-Power  Microwave  Sources  and  Technologies”  [1], 

Attachments: 

X.  Yu,  J.P.  Verboncoeur  and  D.L.  Bruhwiler,  “A  Radiation  Damping  Algorithm  for  Particle 
Simulation,”  unpublished  report. 

J.P.  Verboncoeur  and  K.L.  Cartwright,  “Accuracy  Analysis  of  the  PIC  Method”,  Bull.  Am.  Phy. 
Soc.  45  (2000). 

D.L.  Bruhwiler,  “New  Particle-in-Cell  Algorithms  for  Improving  Accuracy  and  Reducing 
Noise,”  presented  at  the  AFOSR  Electromagnetics  Workshop  (San  Antonio,  January,  2001). 

D.L.  Bruhwiler,  K.  Luetkemeyer,  J.R.  Cary,  D.A.  Dimitrov  and  P.  Stoltz,  “Error  Analysis,  Self¬ 
force  Effects  &  Interactive  Visualization  for  Particle-in-Cell  Codes,”  presented  at  the 
Particle-in-Cell  &  Finite  Difference  Time  Domain  Workshop  (Albuquerque,  August,  2001). 

J.P.  Verboncoeur,  “Symmetric  spline  weighting  for  charge  and  current  density  in  particle 
simulation”,  accepted  for  publication  in  J.  Comp.  Phys.  (2001). 

J.  P.  Verboncoeur,  “Digital  filtering  in  radial  coordinates  for  particle  simulation”,  to  be  submitted 
to  J.  Comp.  Phys.  (2001). 

I.D.  Boyd  and  J.  Fang,  “Final  Technical  Report  -  Fluctuation  Reduction  Algorithms  for  Particle 
Simulations  of  Plasmas.” 

I.D.  Boyd  and  J.  Fang,  “Fluctuation  Reduction  Algorithms  for  Particle  Simulations  of  Plasmas” 
presented  at  the  Particle-in-Cell  &  Finite  Difference  Time  Domain  Workshop  (Albuquerque, 
August,  2001). 

Interactions  /  Transitions: 

The  principal  investigator  Dr.  David  Bruhwiler  and  Prof.  John  Cary  of  the  University  of 
Colorado  and  Tech-X  Corp.  both  attended  the  “AFOSR  Electromagnetics  Workshop,”  organized 
by  Dr.  Arje  Nachman  and  held  January  11-13  in  San  Antonio.  Dr.  Bruhwiler  made  a 
presentation  [10],  which  is  attached  to  this  report. 

The  principal  investigator  Dr.  David  Bruhwiler  and  the  Pi’s  for  the  two  subcontracts.  Prof. 
Ned  Birdsall  and  Prof.  Iain  Boyd,  attended  the  “Particle-in-Cell  /  Finite  Difference  Time 
Domain”  workshop,  organized  by  Dr.  Arje  Nachman  and  Dr.  Bob  Peterkin  and  held  August  21- 
22,  2001  at  Kirtland  AFB.  Dr.  Bruhwiler  made  a  presentation  [11],  which  is  attached  to  this 
report.  Prof.  Boyd  made  a  presentation  [14],  which  is  also  attached  to  this  report. 

The  most  recent  version  of  the  XOOPIC  code  has  been  made  available  to  Dr.  Keith 
Cartwright  of  Kirtland  AFB. 
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Analysis  of  errors  in  the  fundamental  2-D  PIC  algorithm 

-  methodology 
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Error  Analysis  of  the  Generic  PIC  Algorithm 

Work  done  primarily  by  K.  Cartwright  &  J.  Verboncoeur 
Motivation 
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Error  Flow  Schematics 
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Mathematica  is  Used  to  Handle  the  Algebraic  Complexity 
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Charge  Deposition  for  Electrostatic  PIC 


Shape  Function  determined  by  Deposition  Algorithm 
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Cancellation  of  Errors  can  Occur  for  Linear  Deposition 
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Error  Cancellation  not  Possible  for  Electromagnetic  Case 
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Analysis  of  the  Self-Force  Effect  for  Electrons 
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Testing  Self-Force  Implementation  in  XOOPIC 
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Testing  Self-Force  Implementation  in  XOOPIC 

10  GeV  electron  in  a  75  Tesla  magnetic  field 

-  energy  loss  of  single  particles  agrees  well  with  theory 
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Interactive  GUI  for  Particle-in-Cell  Codes 
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Qt  —  a  cross-platform  C++  windowing  toolkit 
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Error  Flow  Schematics 
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Figure  1:  PIC  flow  schematic. 
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Figure  2:  Electrostatic  PIC  error  flow  schematic. 
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Symmetric  current  weightings  are  not  charge  conserving  (zero  order  in  Poisson 
Equation)  NGP  interpolation  in  the  direction  of  the  current  component  and 
linear  interpolation  in  the  transverse  direction  is. 
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Here  a  =mod3  (a  +  1) ,  a  =mod3  {ol  +  2)  are  right-handed  permutations  of  the 
index  a. 
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Here  C  and  C  have  the  units  of  capacitance  and  inductance,  respectively,  and 
store  the  geometric  constants  of  the  mesh  and  the  permittivity  and  permeability. 
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1.  Introduction 

The  Particle  in  Cell  (PIC)  technique  is  a  widely  used  particle  method  for  simulating  plasma  dynamics 
[1].  As  with  many  particle  techniques,  the  solutions  provided  by  the  PIC  method  can  suffer  from  significant 
statistical  fluctuations  under  certain  physical  conditions  of  interest.  In  principle,  these  fluctuations  can  be 
reduced  by  a  combination  of  using  an  increased  number  of  particles  in  a  computation  and  by  time-averaging 
the  solution  over  a  large  number  of  iterations.  Both  of  these  procedures,  however,  lead  to  significant  increases 
in  the  overall  solution  time  and  may  not  work  for  all  conditions. 

The  direct  simulation  Monte  Carlo  (DSMC)  is  a  widely  used  particle  method  for  simulating  neutral, 
rarefied  gas  flows  [2].  The  solutions  generated  by  the  DSMC  technique  also  suffer  from  problems  of  statistical 
noise,  particularly  for  low-speed  flows  that  require  a  large  number  of  samples  to  reduce  the  noise.  Such  a 
simulation  is  extremely  time-consuming  and  beyond  the  capabilities  of  current  computers  [3].  To  solve  this 
problem,  a  new  particle  method  called  the  Information  Preservation  (IP)  algorithm  was  developed  [4-6]. 

In  the  present  report,  we  investigate  whether  the  IP  algorithm  can  be  developed  to  reduce  the  statistical 
noise  in  PIC  simulations.  A  hybrid  PIC/IP  algorithm  is  described  first.  Then  the  hybrid  method  is  assessed  in 
comparison  to  the  standard  PIC  approach  through  its  application  to  three  standard  flows:  (1)  an  unperturbed 
plasma  beam;  (2)  a  perturbed  plasma  beam  (landau  damping);  and  (3)  a  two-stream  instability  problem. 
Finally,  some  conclusions  and  suggestions  for  further  study  are  provided. 

2.  Combination  of  IP  and  PIC  methods 

The  basic  premise  of  the  IP  technique  is  to  associate  both  microscopic  and  macroscopic  information  with 
each  particle  in  a  computer  simulation  [4-6] .  The  usual  PIC  microscopic  steps  are  applied  to  the  microscopic 
properties,  and  simple  conservation  laws  are  applied  to  the  macroscopic  properties.  Macroscopic  quantities 
are  obtained  through  time  or  ensemble  averaging  of  the  information. 

To  clarify  why  the  IP  method  works,  compare  the  formulas  used  by  PIC  and  IP  to  compute  the  mean 
velocity  u  in  a  cell  for  a  warm  plasma.  In  PIC,  the  mean  velocity  is  calculated  from: 


1 


where  N  is  the  sample  size  of  simulated  particles  in  the  cell,  Uk  is  the  velocity  vector  of  particle  k  that 
includes  a  part  ck  due  to  random  thermal  motion.  In  the  IP  method,  the  average  flow  velocity  is  calculated 
using: 

«  =  *£D*  <2> 

k=  1 

where  Uk  is  the  preserved  (macroscopic)  velocity  of  particle  k . 

The  PIC  method  stores  uk,  and  uses  it  to  compute  the  particle  trajectory  and  u  according  to  Eq.  (1). 
The  IP  method  stores  both  uk  and  Uk.  It  uses  the  former  to  compute  the  particle  trajectory,  and  uses  the 
latter  to  compute  u  according  to  Eq.  (2).  Note  that  the  macroscopic  velocity  given  by  the  IP  formula  (2)  is 
the  exact  value  of  u  for  any  sample  size,  whereas  there  is  the  following  explicit  fluctuation  term  in  the  PIC 
formula  (1)  that  decreases  with  the  square  root  of  the  sample  size,  N: 

a=JjY,Ck  (3) 

k~  1 


The  IP  technique  has  proven  very  effective  in  reducing  statistical  scatter  [4-6]  in  particle  simulations 
of  low-speed  rarefied  gas  flows.  For  example,  for  benchmark  problems  such  as  the  Rayleigh  flow  with  a 
characteristic  velocity  of  1  m/s,  the  meaningful  velocity  and  surface  shear  stress  distributions  were  obtained 
using  the  IP  technique  from  a  sample  size  of  103-104.  This  sample  was  four  orders  of  magnitude  smaller 
than  that  required  by  DSMC  (about  108)  [4]  resulting  in  a  tremendous  reduction  in  CPU  time  with  the  IP 
approach. 

An  implementation  of  the  IP  algorithm  for  the  PIC  particle  method  may  be  summarized  as  follows: 

1.  Assign  each  simulated  particle  k  an  information  velocity  Uk,  density  pk,  and  temperature  Tk,  and  a 
microscopic  (PIC)  velocity  uk,  and  assign  each  computational  cell  l  a  macroscopic  velocity  Uj,  density 
pi ,  and  temperature  7). 

2.  Initially,  set  Ui,  pi,  and  of  each  cell  according  to  the  initial  flow  conditions,  and  set  Uk,  pk >  and  2* 
of  all  simulated  particles  in  the  cell  to  be  equal  to  the  cell  values  Ui,  p\,  and  7). 

3.  Move  the  simulated  particles  according  to  the  microscopic  velocities  Uk  using  the  same  algorithms  and 
models  as  the  PIC  method  that  are  described  in  detail  in  Ref.  1. 

4.  In  a  time  step  At,  the  information  velocity  of  a  simulated  particle  may  be  changed  due  to  acceleration 
by  external  forces.  The  sum  of  the  forces  acting  on  cell  l  is 

2(f> 

Fl  =  '£plAA{ei  +  fiAVl  (3) 

j= 1 

where  p[  is  the  gas  pressure  acting  on  surface  j  of  cell  l  with  area  A A\  and  normal  direction  e\,  <j>  is  the 
flow  dimension,  and  //  is  the  volume  force  acting  on  the  cell  with  volume  AV/.  The  term  $  is  computed 
using  the  pressures  in  the  cells  on  either  side  of  the  surface  j .  The  cell  pressures  are  calculated  from  the 
values  of  pi  and  7/  using  the  ideal  gas  equation  of  state.  The  resultant  acceleration  is: 


(4) 


2 


which  contributes  a  velocity  increment  aiAt  to  every  simulated  particle  in  the  cell  during  the  time  step. 

5.  Update  Ui,  pi  and  Tt  each  time  step  as  follows: 

(a)  The  velocity  is  updated  using 

U,  =  wUnew  +  (1  -  w)Uold  (6) 


where  a;  is  a  relaxation  factor  ranging  between  0  and  1,  UQid  is  the  value  of  Uj  from  the  last  time 
step,  Unew  is  the  arithmetic  mean  of  the  information  velocities  of  all  the  simulated  particles  in  the 
cell  during  this  time  step: 


Unew  — 


Eti  mk 

where  m*  is  the  mass  of  particle  k.  For  large  N ,  u  may  be  set  to  1. 

(b)  The  density  may  be  updated  using  the  updated  cell  velocity  and  the  mass  conservation  equation 


(7) 


dp  W) 

dt  dxi 


(8) 


This  equation  is  differenced  using  the  following  implicit  scheme: 


ni+ 1  _  ni  ni+1ui+1  —  flH-lr/H-l 
Pi  Pi  ,  Pl+\Ul+l  Pj-\Ul-\  _  Q 

At  2Ax 

where  superscript  i  denotes  the  values  at  time  step  i. 

(c)  From  the  kinetic  definition  of  temperature,  we  have 


m(c?  -  Ci2) 

J-l,i  ~  - , - 

kjs 


(id 


(12) 


where  a  bar  over  a  quantity  denotes  the  average  value  over  all  the  particles  in  the  sample,  the  subscript  i 
(=1,2,3)  denotes  the  direction,  and  ks  is  the  Boltzmann  constant. 

6.  Update  the  information  densities  of  all  the  simulated  particles  in  cell  l: 


Pk  =  Pi 


(13) 


7.  Compute  macroscopic  quantities  through  statistical  averaging  of  the  information  quantities. 

8.  For  steady  flows,  repeat  steps  3-6  until  the  flow  reaches  a  steady  state.  Then  repeat  steps  (3-7)  to 
sample.  For  unsteady  flows,  repeat  the  steps  (2-7)  until  the  end  of  the  evolution  period. 

The  IP  steps  described  above  have  been  incorporated  into  XES1,  an  electrostatic,  one-dimensional  code 
developed  by  the  Plasma  Theory  and  Simulation  Group  headed  by  Professor  Birdsall  at  UC-Berkeley.  In 
the  following  section,  results  obtained  with  the  new  PIC-IP  formulation  are  compared  with  those  obtained 
with  the  original  PIC  form  of  the  XES1  code. 


3.  Results 

Three  different  test  cases  involving  warm  electrons  are  considered  to  allow  assessment  of  the  usefulness 
and  accuracy  of  the  PIC-IP  scheme:  (1)  an  unperturbed  stream  of  electrons;  (2)  a  perturbed  stream  of 
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electrons  exhibiting  Landau  damping;  and  (3)  a  two-stream  instability  flow.  These  flows  are  used  to  compare 
the  PIC  and  PIC-IP  results  in  terms  of  their  statistical  fluctuations,  their  physical  accuracy,  and  their  overall 
solution  time.  The  codes  use  normalized  input  and  output  data  in  which  time  is  normalized  by  the  plasma 
frequency  and  distance  is  normalized  by  the  Debye  length. 

3.1  Unperturbed  Stream  of  Electrons 

The  problem  is  initialized  with  a  warm  stream  of  electrons  as  follows:  stream  velocity,  Vi=3.0;  length  of 
system,  l=2ir;  there  is  no  perturbation  applied  to  either  the  particle  velocities  or  positions,  and  there  is  no 
applied  electric  field.  The  numerical  parameters  are:  time  step,  At=0.2;  number  of  computational  cells=32; 
number  of  iterations=10,000;  number  of  particles  per  cell  is  varied. 

For  this  unperturbed  system,  the  steady  state  solution  is  a  flat  profile  at  the  initial  velocity,  vx>  In 
Fig.  la,  PIC  profiles  are  shown  at  the  beginning  and  the  end  of  the  simulation  that  are  obtained  using 
50  particles  per  cell.  These  results  indicate  that  the  simulation  does  provide  the  expected  steady  state 
result,  but  the  data  are  accompanied  by  a  significant  degree  of  statistical  fluctuation.  These  fluctuations  are 
quantified  in  terms  of  the  standard  deviation  and  listed  in  Table  1  for  a  number  of  different  simulations.  The 
effect  of  increasing  the  number  of  particles  per  cell  can  be  seen  in  Fig.  lb.  The  fluctuations  are  successively 
reduced  as  the  number  of  particles  is  increased  such  that  a  one  order  of  magnitude  reduction  in  standard 
deviation  is  achieved  by  increasing  the  total  number  of  particles  by  two  orders  of  magnitude,  as  expected 
from  Eq.  (3). 

A  technique  commonly  employed  in  PIC  simulations  to  reduce  statistical  fluctuations  is  the  so-called 
“quiet  start”  in  which  the  initial  values  of  the  particle  positions  are  randomized  according  to  the  initial 
particle  velocities  so  as  to  more  uniformly  distribute  the  particles  in  velocity-position  space.  The  effects  of 
using  a  quiet  start  on  the  unperturbed  stream  are  shown  in  Fig.  lc.  From  the  figure  and  from  Table  1  we 
can  see  that  while  the  quiet  start  does  initially  reduce  the  level  of  statistical  fluctuations  at  the  start  of  the 
simulation,  the  fluctuations  levels  ultimately  reach  about  the  same  level  as  those  obtained  with  the  same 
number  of  particles  without  the  quiet  start. 

Finally,  in  Fig.  Id,  the  PIC-IP  solutions  at  the  beginning  and  end  of  the  simulation  are  shown  and 
compared  with  the  corresponding  PIC  results.  As  confirmed  in  Table  1,  the  level  of  fluctuations  achieved 
at  the  end  of  the  PIC-IP  simulation  using  500  particles  per  cell  is  smaller  than  that  obtained  with  the  PIC 
method  using  5,000  particles  per  cell.  For  completeness,  the  PIC-IP  and  PIC  results  are  also  compared  for 
the  electric  field,  and  charge  density  p  in  Figs,  le  and  If,  respectively.  These  show  essentially  the  same 
trends  as  the  velocity  profiles. 

Included  in  Table  1  are  the  CPU  times  for  all  the  simulations.  These  are  the  overall  run  times  including 
all  data  input  and  output.  The  data  show  that  the  PIC-IP  method  using  500  particles  per  cell  obtains 
solutions  with  smaller  statistical  fluctuations  than  the  PIC  method  with  5,000  particles  per  cell  in  a  CPU 
time  that  is  smaller  by  a  factor  of  more  than  eight. 

Thus,  the  simple  test  case  of  an  unperturbed  stream  shows  several  important  results.  PIC  simulations 
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are  statistically  noisy,  and  these  fluctuations  can  only  be  suppressed  by  employing  a  significantly  larger 
number  of  particles.  The  PIC-IP  method  successfully  reduces  these  fluctuations  and  achieves  the  desired 
physical  solution,  allowing  significant  reductions  in  computational  cost. 

3.2  Perturbed  Stream  of  Electrons  (Landau  Damping) 

The  flow  of  the  previous  section  is  modified  by  perturbing  the  initial  particle  positions  according  to  xj 
cos(27rx/Z)  where  xi=0.1,  and  changing  the  mean  velocity  to  zero. 

Profiles  at  three  different  times  of  charge  density,  electric  field,  and  stream  velocity  computed  using 
both  the  PIC  and  PIC-IP  algorithms  are  shown  in  Figs.  2a-2c,  respectively.  In  each  simulation,  an  average 
of  500  particles  per  cell  is  used.  In  Fig.  2a,  the  density  profile  shows  a  gradual  damping  from  the  initial 
perturbation,  and  this  is  also  seen  in  the  electric  field  profile.  The  stream  velocity  profile  shows  an  interesting 
feature  where  a  sinusoidal  profile  is  generated  after  10  iterations  that  eventually  is  damped  back  to  a  nearly 
flat  profile.  In  all  of  these  plots,  it  can  be  seen  that  there  is  significant  fluctuations  in  the  PIC  results,  and 
that  the  PIC  data  tend  to  fluctuate  about  the  smoother  PIC-IP  profiles.  In  Fig.  2d,  the  velocity  profile 
obtained  using  5,000  particles  per  cell  with  PIC  is  compared  with  the  PIC-IP  solution  obtained  with  500 
particles  per  cell.  In  comparison  with  the  PIC  data  shown  in  Fig.  2c,  the  PIC  data  shown  in  Fig.  2d  exhibit 
smaller  fluctuations  and  more  clearly  fluctuate  about  the  PIC-IP  solution. 

In  Fig.  2e,  the  energy  field  histories  are  shown  for  the  PIC  and  PIC-IP  simulations  employing  500 
particles  per  cell.  The  PIC  results  show  the  expected  behavior  of  the  energy  being  reduced  by  more  than 
an  order  of  magnitude  before  undergoing  regular  oscillations.  As  described  in  Ref.  1,  further  reduction  in 
the  field  energy  can  be  achieved  with  the  PIC  technique  by  dividing  the  velocity  distribution  function  into 
different  species.  The  values  of  the  field  energy  predicted  by  the  PIC-IP  algorithm  are  in  good  general 
agreement  with  the  PIC  results  although  they  appear  to  level  out  at  a  slightly  lower  energy  and  oscillate 
with  larger  amplitude.  It  should  be  noted  that  these  computations  are  performed  without  using  the  quiet 
start  feature. 

These  comparisons  appear  to  indicate  that  the  PIC-IP  algorithm  is  again  able  to  significantly  reduce 
the  statistical  fluctuations  of  the  PIC  method  while  maintaining  physical  accuracy.  The  overall  run  time 
for  the  PIC-IP  simulation  using  500  particles  per  cell  is  about  a  factor  of  seven  smaller  than  that  for  the 
PIC  simulation  employing  5,000  particles  per  cell.  The  reduction  in  computational  time  obtained  using  the 
PIC-IP  approach  is  slightly  smaller  than  that  for  the  simple  stream  considered  in  the  previous  section  only 
because  the  number  of  iterations  is  smaller  here  and  this  leads  to  a  proportionally  larger  time  spent  in  the 
input  and  output  of  data. 

3.3  Two-Stream  Instability 

Consider  two  opposing  streams  of  electrons  with  equal  number  density,  n0,  equal  temperature,  T,  and 
equal  velocities  moving  in  opposite  directions,  ±«i.  There  is  a  background  of  immobile  positive  charges 
with  a  total  density  of  2no-  When  the  two  streams  traverse  one  wavelength  during  one  plasma  oscillation, 
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a  density  perturbation  in  one  of  the  streams  is  increased  by  the  electrostatic  forces  created  by  bunching  of 
particles  in  the  other  stream,  and  vice  versa.  Hence,  the  perturbation  is  observed  to  grow  exponentially  with 
time,  and  the  system  is  unstable.  This  type  of  system  is  studied  here  using  the  following  physical  parameters: 
stream  velocity,  dhl.O;  thermal  velocity  (quiet  start),  u*=0.5;  length  of  system,  I=107r;  initial  particle 
positions  are  perturbed  according  to  xi  cos(2irx/l)  where  x\  =±0.001.  The  numerical  parameters  are:  a  time 
step,  At=0.1;  number  of  computational  cells— 128;  number  of  iterations=300;  initial  number  of  particles  per 
cell  is  varied  between  64  and  4096. 

In  Fig.  3a,  the  charge  density  profile  after  completion  of  one  iteration  of  the  PIC  solution  is  shown  for 
various  values  of  average  number  of  particles  per  cell.  Clearly,  the  initial  charge  density  profile  is  strongly 
affected  by  the  number  of  particles  employed  with  significant  reductions  in  statistical  fluctuations  obtained 
when  a  larger  numbers  of  particles  is  employed.  Similar  results  are  found  for  the  electric  field  and  the  mean 
velocity,  as  shown  in  Figs.  3b  and  3c,  respectively.  In  the  case  of  the  stream  velocity,  only  the  profile  for  the 
forward  moving  population  of  electrons  is  shown.  The  population  of  backward  moving  electrons  has  identical 
behavior.  In  Fig.  3d,  the  PIC  (4096  particles  per  cell)  and  PIC-IP  (512  particles  per  cell)  profiles  of  charge 
density  after  one  iteration  are  compared.  Clearly,  the  PIC-IP  algorithm  is  again  capable  of  significantly 
reducing  the  fluctuations. 

In  Figs.  3e  through  3g,  the  density,  electric  field,  and  mean  velocity  profiles  obtained  at  the  end  of 
the  simulation  using  various  numbers  of  particles  per  cell  with  the  PIC  method  are  shown.  It  is  disturbing 
that  the  final  solutions  for  all  variables  depend  so  strongly  on  the  number  of  particles  employed.  These 
same  profiles  are  compared  for  the  PIC  solution  employing  4096  particles  per  cell  and  the  PIC-IP  solution 
employing  512  particles  per  cell  in  Figs.  3h  through  3j.  While  there  are  differences  between  the  solutions, 
they  are  also  of  a  more  similar  nature  than  the  various  PIC  solutions  shown  in  Figs.  3e-3g.  While  these 
types  of  flow  require  more  careful  study,  the  PIC-IP  results  are  encouraging  in  that  solutions  similar  to  the 
full  PIC  method  can  be  obtained  with  a  factor  of  eight  fewer  particles  (corresponding  to  an  overall  reduction 
in  computation  time  of  a  factor  of  six). 

4.  Conclusions 

In  this  report,  we  studied  the  possibility  of  applying  the  IP  technique  to  the  PIC  method.  A  hybrid 
PIC/IP  algorithm  was  proposed,  and  was  applied  to  simulate  three  different  one-dimensional,  electro-static 
problems:  (1)  an  unperturbed  electron  stream;  (b)  a  perturbed  electron  stream  (landau  damping);  and  (3) 
a  two-stream  plasma  instability.  In  each  case,  the  same  general  conclusions  were  reached.  First,  the  PIC 
solutions  exhibited  strong  statistical  fluctuations  when  small  numbers  of  particles  were  employed.  Second, 
PIC-IP  solutions  employing  factors  of  8  to  10  fewer  particles  exhibited  reduced  statistical  scatter  and  good 
general  agreement  with  the  best  PIC  solutions. 

This  was  by  no  means  a  complete  study  and  there  are  several  potential  directions  for  further  work: 

(1)  A  careful  study  should  be  made  of  the  implementation  and  effect  in  the  PIC-IP  algorithm  of  using 

additional  fluctuation  reduction  techniques  such  as  the  quiet  start  commonly  used  in  PIC. 
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(2)  An  electron  energy  equation  should  be  solved  to  compute  changes  in  the  macroscopic  electron  temper¬ 
ature. 

(3)  Cold  plasma  conditions  should  be  studied. 

(4)  The  PIC-IP  algorithm  should  be  extended  to  include  magnetic  fields,  and  multi-dimensional  flows. 
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Table  1.  Numerical  properties  for  simulations  of  an  unperturbed  electron  beam. 
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Fig.  le.  Electric  field  profile  for  unperturbed  stream. 
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Fig.  If.  Density  profile  for  unperturbed  stream. 
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Fig.  2b.  Electric  field  profile  for  perturbed  stream. 
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Fig.  2e.  Field  energy  history  for  perturbed  stream. 


Fig.  3a.  Initial  density  profiles  for  warm,  two-stream  instability. 
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Fig.  3d.  Initial  density  profiles  for  warm,  two-stream  instability. 


Fig.  3e.  Final  density  profiles  for  warm,  two-stream  instability. 
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Abstract 

Radiation  damping  is  a  phenomenon  that  occurs  when  a  charged 
particle  undergoing  a  large  acceleration  emits  radiation.  While  this 
occurs  in  a  linear  system  only  with  extreme  accelerations,  the  effect 
is  more  pronounced  for  a  magnetized  acceleration  where  the  radial 
acceleration  induced  by  the  Lorentz  force  is  sufficient  to  generate  syn¬ 
chrotron  radiation.  The  effect  is  only  measurable  in  highly  relativistic 
charged  particles.  The  radiation  damping  results  in  loss  of  energy 
and  momentum,  so  that  an  electron  orbit  will  decay  in  time.  In  this 
work,  we  consider  a  relativistic  electron  orbiting  in  a  uniform  magnetic 
field.  A  model  following  Jackson  [1]  is  implemented  in  the  oarticle- 
in-cell  code  XOOPIC  [2].  The  model  includes  the  momentum  and 
energy  loss  effects,  but  neglects  the  propagation  of  the  emitted  radi¬ 
ation.  Simulation  results  are  obtained  for  a  10  GeV  electron  in  a  75 
T  field,  and  the  temporal  decay  of  the  energy  compares  favorably  to 
theory. 

*TechX  Corporation,  1280  28th  St.,  Suite  2,  Boulder,  CO  80303 
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1  Introduction 


Radiation  damping  is  a  physical  phenomenon  in  which  a  charged  particle 
undergoing  rapid  acceleration  emits  radiation,  resulting  in  loss  of  particle 
energy  and  momentum.  The  rapid  acceleration  may  be  related  to  forces  from 
electromagnetic  fields  or  from  collisions.  This  effect,  sometimes  referred  to 
as  self-force,  is  only  of  significance  for  highly  relativistic  charged  particles 
undergoing  extreme  acceleration. 

While  collisions  can  also  result  in  significant  accelerations,  the  energy  loss 
associated  with  a  collision  can  be  treated  as  part  of  the  collision  mechanics; 
we  are  seeking  effects  here  that  modify  the  traditional  Newton-Lorentz  equa¬ 
tions  of  motion.  In  addition  to  collisions,  there  are  two  general  classical  cases 
where  radiation  damping  is  exhibited.  One  is  linear  acceleration  of  electrons 
in  a  strong  electric  field.  The  other  is  synchrotron  radiation  emitted  by  elec¬ 
trons  accelerating  in  due  to  circular  motion  about  a  magnetic  field.  For  the 
linear  acceleration  case,  the  energy  loss  is  small  for  most  achievable  param¬ 
eters,  and  is  thus  not  directly  addressed  here.  Synchrotron  radiation,  on  the 
other  hand,  is  readily  observed  and  is  a  significant  factor  in  the  design  of 
particle  accelerators.  In  this  research,  the  case  of  a  particle  accelerating  in 
circular  motion  about  a  uniform  static  magnetic  field  is  considered,  since  this 
creates  a  classical  acceleration  of  magnitude  sufficient  to  induce  significant 
radiation  damping. 

A  macroscopic  view  of  radiation  damping  is  discussed  by  Jackson  [1], 
including  the  linear  acceleration  and  synchrotron  cases  as  well  as  the  angle 
and  spectrum  of  the  emitted  radiation.  Yaghjian  [3]  presents  a  detailed  study 
for  a  spherically  shaped  particle,  and  Rossi  [2]  presents  a  brief  discussion 
on  the  topic.  The  general  conclusions  of  the  authors  are  similar,  although 
the  complexity  and  level  of  detail  varies  greatly.  This  study  focuses  on  the 
Jackson  model  due  to  its  simplicity  and  its  applicability  to  the  PIC  model. 
The  authors  are  not  aware  of  any  other  computational  work  in  the  literature 
attempting  to  address  radiation  damping. 

The  objective  of  the  research  was  to  implement  and  analyze  a  compu¬ 
tationally  efficient  correction  to  the  standard  particle-in-cell  (PIC)  scheme 
to  account  for  the  particle  kinetics  when  radiation  damping  becomes  non- 
negligible.  This  work  makes  the  PIC  algorithms  more  accurate  for  modeling 
high-energy  particles  undergoing  rapid  accelerations  by  modifying  the  equa¬ 
tions  of  motion  to  account  for  the  loss  of  energy.  The  radiated  energy  is  not 
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treated  here;  in  particular,  the  spectral  and  directional  distribution  of  the 
radiated  energy  is  not  considered,  and  the  interaction  of  the  radiated  fields 
with  the  particle  is  also  neglected;  the  radiated  energy  is  simply  treated  as 
lost  energy. 

The  theoretical  model  employed  is  discussed  in  Section  2,  including  a 
discussion  of  the  computational  resources  required  to  include  the  radiated 
energy  and  fully  resolve  the  generated  fields  and  their  effects  on  the  particle. 


2  Radiation  Damping  Model 
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Figure  1:  The  Model  used  for  the  algorithm. 


In  this  section,  a  we  start  with  a  general  model  for  radiation  damping 
from  Jackson  [1],  The  linear  acceleration  is  discussed,  and  it  is  demonstrated 
that  the  effect  is  negligible  for  this  case.  Next,  an  electron  accelerating  in 
circular  motion  in  a  static  homogenous  magnetic  field  is  considered,  again 
following  Jackson.  The  power  loss  is  described  for  the  synchrotron  case,  and 
the  spectral  and  direction  emission  of  radiation  is  discussed,  as  well  as  the 
ramifications  of  modeling  the  emitted  radiation  in  a  PIC  code.  Finally,  the 
radiation  damping  algorithm  is  discretized  for  incorporation  into  the  PIC 
model. 

For  an  electron  undergoing  a  rapid  acceleration,  Jackson  related  the  ra¬ 
diated  power  to  the  derivative  of  the  electron  momentum: 


P  =  - 


3  47 TEom^c3 


(1) 
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where  e  is  the  electronic  charge,  m  is  the  electron  rest  mass,  c  is  the  speed 
of  light,  £o  is  the  permittivity  of  free  space,  dr  =  dt/ 7  is  the  proper  time 
element,  7  =  (1  —  p)~1^2 ,  p  =  7 mv,  /3  =  v/c,  and  v  is  the  electron  velocity. 
For  linear  acceleration,  the  corresponding  result  is  written: 


P  =  2  el  f— V 

3 47r£0ra2c3  \dt ) 


(2) 


It  is  instructive  to  note  that  in  the  linear  acceleration  case,  Jackson  showed 
that  the  radiation  losses  are  negligible  compared  to  the  effects  of  the  external 
forces  acting  on  the  electron  when  the  change  in  energy  per  unit  distance 
satisfies 

dE 
dx 


( 2  e2  y 
<  \  3  47T£o m2c4/?  i 


(3) 


which  becomes  2.7  x  1020  V/m  for  highly  relativistic  electrons.  Since  this 
is  far  above  the  gradients  of  107  V/m  present  in  linear  accerators,  radiation 
damping  is  a  negligible  effect  for  linear  devices.  We  therefore  consider  only 
the  synchrotron  model  in  this  research. 

Consider  an  electron  in  circular  motion  in  a  static  magnetic  field  as  shown 
in  Fig.  1.  For  synchrotron  motion,  the  radiated  power  becomes  [1]: 


2  e2 

3  47T£0 t\ 


4„,4 


Pi 


(4) 


and  the  relativistic  Larmor  radius  is  given  by: 


v  _  7 m/3c 
ujc  eB 


(5) 


where  B  is  the  external  static  magnetic  field. 

The  spectral  content  of  the  radiation  power  is  a  complicated  function. 
The  radiation  is  strongly  polarized  in  the  plane  of  motion,  with  increasing 
confinement  to  the  plane  of  motion  with  increasing  frequency.  The  radiation 
power  is  confined  to  increasingly  smaller  angles  at  higher  frequencies,  as  well 
as  with  increasing  particle  energy.  Jackson  derived  the  energy  radiated  per 
unit  frequency  per  unit  solid  angle  for  the  synchrotron  case: 


d2/ 

diodQ 


(*f/3  R)  +  R>) 


(0 
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where  6  is  the  angle  from  the  plane  of  circular  motion,  K  is  the  modified 
Bessel  function,  and  f  is  given  by: 

(7) 

The  radiation  is  negligible  beyond  the  critical  frequency, 

Merit  =  373— .  (8) 

tl 

The  number  of  harmonics  present  in  the  radiation  spectrum  for  periodic 
motion  of  angular  frequency  lo  is  given  by  nCTit  =  37s.  The  critical  angle 
beyond  which  the  radiated  power  drops  to  1/e  of  the  peak  value  is  given  by 

9crit  ~  7/^ 

From  Eq.  9,  it  is  evident  that  the  radiated  power  is  confined  to  increasingly 
smaller  angles  at  higher  frequencies,  as  well  as  increasing  particle  energy. 
The  radiation  is  emitted  from  the  particle  undergoing  circular  motion  in  a 
pattern  similar  to  the  sweep  of  a  lighthouse  beam,  with  the  focus  in  the 
direction  of  the  velocity  vector. 

To  assess  the  time  and  space  scales  required  to  resolve  the  emitted  ra¬ 
diation,  consider  a  1  GeV  electron  orbiting  a  10  T  static  magnetic  field. 
The  relativistic  momentum  factor  is  7  =  1952,  and  the  velocity  is  effectively 
v  =  c.  The  gyroradius  is  77  =  0.334  m.  The  power  radiated  is  P  =  6  x  10-6 
W,  so  the  electron  loses  about  1%  of  its  energy  in  270  ns.  The  critical  fre¬ 
quency  Mcrn / 2tt  =  3.2  x  1018  Hz,  with  a  corresponding  critical  wavelength  of 
A crit  =  9-4  x  10-11  m.  In  order  to  resolve  the  wavelength,  the  mesh  size  re¬ 
quirement  becomes  Ax  <  A/4  =  2.4  x  10-11  m.  For  a  nominal  system  length 
of  L  =  3 77  in  each  dimension,  the  required  number  of  cells  in  each  dimension 
is  4.2  x  1010,  or  1.7  x  1021  cells  in  a  two-dimensional  simulation.  The  storage 
required  for  this  spatial  resolution  exceeds  the  memory  capacity  of  exisiting 
computers.  Furthermore,  the  timestep  required  to  satisfy  the  Courant  con¬ 
dition  is  At  <  Ax/y/2 c  =  5.7  x  10~20  s.  The  cyclotron  period  is  7  ns,  so  one 
orbit  would  require  1.2  x  1011  timesteps,  and  running  to  the  1%  energy  loss 
would  require  4.7  x  1012  timesteps.  Resolving  the  critical  frequency  requires 
a  timestep  At  <  3.1  x  10~19  s,  so  frequency  resolution  is  even  more  unten¬ 
able.  It  is  clear  from  these  space  and  time  constraints  that  it  will  remain 
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impossible  for  the  forseeable  future  to  simulate  the  self-consistent  emission 
of  radiation  from  even  a  single  particle  within  a  standard  PIC  scheme.  It  is 
equally  impossible  to  simulate  the  emitted  fields  over  a  frequency  spectrum 
sufficient  to  include  a  significant  fraction  of  the  radiated  power. 

The  total  power  radiated  for  the  synchrotron  case  can  be  written  by 
modification  of  Eq.  1: 


P  = 


e2q2 

67T£oC3 


(10) 


where  u  =  yv  and  u  —  \/u  •  u. 

Using  —  P  =  ( dE/dt)rad  for  the  radiation  energy  loss  rate,  where  the 
particle  energy  is  given  by  E  =  yrac2,  we  can  rewrite  the  power  loss: 


± A  1  p_  eV 

dt  J  ,  me 2  6ne0mc5 


(11) 


Writing  the  magnitude  of  the  normalized  momentum  in  terms  of  7: 

u  =  csff  -  1,  (12) 

and  taking  the  time  derivative  to  obtain 

du  =  yc  dy 

dt  y/r^ldt'  {  } 

Combining  Eqs.  11  and  13,  and  applying  the  result  to  radiation  damping: 

2 

/  01  \  /  rim  \  I 

(14) 


du 
dt  , 


7° 


rad 


67T£0  rnc4  ^y2  —  1 


du' 

dt  , 


2  /  ,  X  2 


u  \  ( du\ 


Jcj 


dt  J 


Discretization  goes  here  [David  Bruhwiler] 


3  Analysis  and  Verification  of  the  Algorithm 

In  order  to  verify  that  the  damping  rate  obtained  in  the  simulation  is  correct, 
an  analytical  formula  is  derived  for  the  energy  as  a  function  of  time,  and  is 
compared  to  the  simulation  result  for  the  synchrotron  radiation  case. 
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Figure  2:  Synchrotron  decay  of  single  electron  orbit. 


We  can  write  variables  in  terms  of  energy,  7  =  E/mc 2,  as  well  as  (32  = 
1  —  m2c4/E2.  The  particle  velocity  can  be  written 


v  —  c\  1  — 


so  the  normalized  momentum  becomes 


m2c4 

~E2~} 


u 


I  E2 
m2c 4 


1. 


For  highly  relativistic  particles,  Eq.  16  reduces  to  u  ~  E/mc. 
Then  we  can  write 


du 

dt 


1 


dE 


du  dE 

dE  dt  mc^Jl  —  m2^/^2  ^ 


(15) 


(16) 


(17) 


When  the  energy  is  large  compared  to  the  rest  energy, 

du  1  dE 
_  _ 

dt  me  dt 

Consider  the  equation  of  motion  for  an  electron  in  a  magnetic  field: 

dp  mdu 


dt 


dt 


-ev  x  B. 


(18) 


(19) 
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Here  we  consider  the  synchrotron  case  shown  in  2,  where  the  magnetic  field 
is  perpendicular  to  the  plane  of  motion.  Using  Eq.  15, 


du  eB  m2c4 

dt  m  V  E 2 


In  the  high-energy  limit, 


(20) 


du  eBc 
dt  m 


(21) 


Using  Eq.  21  in  Eq.  10,  we  obtain: 


dE  _ 

1  e2  E2 

’  e2B2 

1 

(  dE\2' 

dt 

67TS0  c3  m2c4 

EoF  om2 

m2c2 

\dt  ) 

To  approximate  the  equation,  compare  with 


1  (dE\  2 
l  dt  ) 


stants  and  parameters  for  electrons,  we  get 

e2U2  =  (1.6  x  10-19)2  x  752  =  2g  x  17 

e0//o  8.85  x  lO”12  x  4?r  x  10~7 


(22) 

Insert  con- 


1  ,dEs2 


1 

(3  x  108)2 


(10-1)2  =  1.1  x  10~19 


where  we  used  approximation  ^  ~  0.1,  which  is  valid  as  seen  from  the 
next  sections.  And  this  approximation  will  keep  being  valid  since  the  energy 
radiation  rate  will  decrease.  So  the  second  term  in  equation(12)  could  be 
dropped  for  approximation. 

Thus,  the  approximated  form  becomes: 


dE  _ 1  e2  E2  e2H2  _  __eW_ 

dt  67re0  c3  m2c4  eofiom?  67reoPo^4c7 


dE  (1.6022  x  10~19)4  x  752  x  E2 

II  ~  _6tt(8.85  x  10-12)24tt  x  10"7  x  (9.1095  x  lO"31)4  x  (2.9979  x  108)7 

HF, 

—  =  -1.3191  x  1016£2  (23) 

dt  v  ’ 

Solving  for, 


8 


(24) 


or, 


E  =  ~ 


1 


E0 


+  1.3191  x  1016t 


4  =  -7-  +  1.3191  x  1016f 

E  E0 


(25) 


4  Results 

As  stated  before,  synchrotron  radiation  could  give  a  good  observation  of  the 
damping.  An  input  file  was  set  up  to  simulate  it.  In  this  simulation,  one 
electron  is  injected  into  a  strong  magnetic  field  with  a  large  initial  energy. 

Attachment  1  gives  the  input  file  used.  And  figure  3  shows  the  x-y  phase 
space  of  the  movement  of  the  electron.  And  the  energy  damping  could  be 
easily  observed  from  the  decreasing  of  its  gyro-radius. 


x-y  phase  space 


Figure  3:  Electron  position  trace  from  simulation. 
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And  the  verification  of  the  simulation  was  done  by  comparing  the  kinetic 
energy  diagnostic  with  the  energy  damping  function  developed  above.  And 
a  comparason  of  the  analytical  energy  function  with  the  energy  loss  function 
from  Jackson’s  was  first  done  to  illustrate  its  correctness. 

4.1  Check  the  energy  formula  against  Jackson’s  one  circle  energy 
dissipation  equation. 

My  formula  (13)  is  in  MKSA.  Check  it  against  Jackson’s  equation  (14.33) 
for  high  energy  electrons  (/?  «  1): 

6E(Mev )  =  8.85  x 

p(meters) 

For  the  first  circulation,  gyration  radius(in  relativistic  form)  is  given  by: 

_  7 mv0  _  1.9570  x  104  x  9.1095  x  lO^31  x  2.9979  x  108 
P~  eB  ~  1.6022  x  10~19  x  75 


0.4448  (m) 


Where  Vo  is  approximated  to  be  light  velocity,  since  VuHft  —  1010eu  as 
given  in  the  input  file.  And  7  is  calculated  as  following: 

Kinetic  energy: 

T  =  (7  —  l)mc2  =  eV 
eV  „  1.6022  x  10“19  x  1010 

7=1  + 


=  1  + 


mr 


9.1095  x  10~31  x  (2.9979  x  108)2 


=  1.9570  x  104 


by 


Then  the  time  for  the  particle  to  move  through  the  first  cycle  is  calculated 

2irp  2irp  2  x  3.14  x  0.4448  nnnnA  .._  q  ,  , 

t  =  —  «  —  = - 5 —  =  9.3224  x  10~9  (s) 

v  c  2.9979  x  108  v  ’ 

Calculate  8 E  by  my  formula: 


~  =  -^+1.3191  xl016t  =  —  1-  -• 
E  E0  1.6  x  10-9 


1. 3191  xl016x  9. 3224x10 


-9 


747113583 


5E  =  E  -  Eo  = 


1 


747113583 
Calculate  8E  by  Jackson’s: 


-  1.6022  x  10"9  =  -2.6152  x  10“10  ( Joules ) 
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SE„ 


-8.85  x  10-2^ 
P 


-8.85  x  10“2  x 
— 1989. 6Mev 


(10)4 

0.4448 


—3.1874  x  10  10  ( jourls ) 


£74 

8 Emir,  =  -8.85  x  10"2 — 


-8.85  x  10-2  x 


(10  -  1.9896)4 
0.4448 


-819.21Meu 

-1.3125  x  10~10  (joules) 


The  result  from  analytical  formula  falls  into  the  interval  between  the 
maximum  and  minimum  energy  change.  Noticing  the  Jackson’s  euqation  is 
also  an  approximation,  this  suggests  the  correctness  of  my  formula. 

4.2  Check  simulation  against  analytical  energy-time  function. 

Pick  two  points  from  KE  diagonistic. 


U  =  1.345  x  10  8 ,  E\  =  0.0001247  x  10~5 
Using  E(t)  function,  E  =  1.2476  x  10“9.  And  the  difference  is  —0.048%. 

t2  =  1.4948  x  10~6,£2  =  4.9986  x  10_u 

Using  E(t)  function,  E  =  4.9214  x  10~H.  And  the  difference  is  0.091%. 
We  are  using  as  defference.  And  all  the  values  above  are  in  MKSA. 
As  we  can  see, the  simulation  result  and  the  analytical  result  are  very  close 
to  each  other.  Attatchment  2,  3,  4,  5  show  the  curves  from  simulation,  from 
analytical  results,  a  combination  of  them  and  the  difference  between  them. 
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5  Conclusion 


As  shown  by  the  analysis  above,  the  radiation  damping  algorithm  imple¬ 
mented  to  PIC  code  is  an  accurate  algorithm  to  approach  the  problem.  This 
makes  PIC  a  more  powerful  tool  in  various  research  work  related  to  plasma 
simulation. 
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Figure  4:  XOOPIC  input  file  for  simulation. 
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1  Secondary  Emission 

The  process  of  electron  impact  secondary  emission  is  a  key  element  of  a  number 
of  processes  related  to  high  power  microwave  tubes.  Secondary  emission  plays 
a  pivotal  role  in  depressed  collectors,  single-  and  two-surface  multipactors,  and 
beam  interception.  In  this  section,  a  secondary  emission  model  is  outlined,  and 
the  implementation  in  the  XOOPIC  PIC-MCC  code  [1]  is  described. 

Electron  impact  secondary  emission  occurs  when  an  electron  impacts  a  sur¬ 
face,  which  may  be  a  conductor  or  a  dielectric,  and  ejects  electrons  from  the 
surface.  In  PIC  codes,  it  is  not  possible  to  model  the  quantum  mechanical 
details  of  the  process,  which  involves  interaction  of  the  incident  electron  with 
conduction  or  valence  band  electrons  in  the  surface  medium,  due  to  the  time  and 
space  scales  involved.  Instead,  it  is  more  efficient  to  employ  a  phenomenological 
model. 


1.1  Energy  and  Angular  Dependence 

This  work  is  based  on  the  secondary  model  due  to  Vaughan  [2]  and  later  ex¬ 
perimentally  verified  by  Shih  [4],  and  improved  by  Vaughan  [3]  and  later  [5]. 
The  secondary  electron  coefficient,  defined  as  the  ratio  of  ejected  to  incident 
electrons,  has  both  energy  and  angular  dependence: 

8  {6,6)  ~  5max0  (wexp(l  —  w))k .  (1) 

Here,  the  incident  energy  is  given  by  6 ,  the  angle  with  respect  to  the  surface 
normal  is  6 ,  kss  is  a  surface  smoothness  parameter  described  below,  A;  is  a 
curve-fit  parameter  also  described  below,  and  £max0  is  the  peak  coefficient,  which 
occurs  at  normal  incidence  at  the  energy  £maxo.  The  energy  dependence  appears 
implicitly  in  the  right  hand  side  of  Eq.  1  through  the  normalized  energy,  w ,  given 
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Figure  1:  Normalized  secondary  emission  coefficient  as  function  of  normalized 
energy  at  normal  incidence,  with  £maxo/£o  —  40,  and  ksw  —  1. 


by: 


_ £-£o _ 

^maxO  (l  4"  k  sw  02/2i r)  -  ’ 


(2) 


where  £q  is  the  secondary  emission  threshold  and  ksw  is  a  surface-smoothness 
parameter  similar  to  kss.  Both  ks$  and  ksw  vary  between  0  for  very  rough 
surfaces  and  2  for  polished  surfaces.  Typical  values  are  close  to  1.  The  exponent 
k  in  Eq.  1  is  given  by: 


f  0.62,  w  <  1 
\  0.25,  w  >  1 


(3) 


The  energy  dependence  of  the  secondary  emission  coefficient  is  shown  in 
Fig.  1.  The  angular  dependence  of  the  secondary  emission  coefficient  is  shown 
in  Fig.  2.  Note  that  the  secondary  emission  coefficient,  which  can  be  obtained 
by  multiplying  the  energy  dependent  part  by  the  angular  dependent  part,  is 
largest  at  £  =  £maxo  and  6  =  90.  The  peak  normal  secondary  emission  in 
copper  is  5max0  =  1.2,  with  £max0  =  400  eV,  and  £q  =  15. 


2 


Figure  3:  Schematic  diagram  of  the  secondary  emission  spectrum  versus  the  ra¬ 
tio  of  the  emitted  energy  of  the  secondary  to  the  incident  energy  of  the  primary. 

1.2  Secondary  Emission  Spectrum 

The  energy  and  angular  distribution  of  the  emitted  secondaries  were  treated  by 
Spangenberg  [6],  and  are  also  treated  by  Vaughan  [2].  The  emission  spectrum 
has  three  regimes,  as  shown  in  Fig.3.  Incident  electrons  reflected  at  the  surface 
of  the  secondary  emitting  material  are  called  reflected  primaries.  Reflected 
primaries  comprise  about  3%  of  the  emitted  electron  population.  The  energy  of 
the  reflected  primary  is  approximately  the  same  as  the  energy  of  the  incident 
primary,  £r  —  £{.  The  primary  is  reflected  specularly  in  angle,  9r  —  — 0i,  where 
9  is  measured  from  the  surface  normal  at  the  point  of  impact. 

Backscattered  primaries  are  electrons  that  impact  the  surface,  and  scatter 
off  of  several  lattice  atoms  and/or  impurities. Typically  backscattered  primaries 
comprise  about  7%  of  the  emitted  electron  population.  These  electrons  re- 
emerge  from  the  impact  surface  with  energies  in  the  range  0  <  £5  <  £.  Within 
this  energy  range,  all  energies  are  taken  as  equally  probable,  £5  =  R£{,  where 
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0  <  R  <  1  is  a  uniformly  distributed  random  number.  The  angle  of  emission  is 
taken  to  be  specular,  just  as  in  the  case  of  the  reflected  primaries.  A  more  de¬ 
tailed  treatment  might  consider  a  distribution  of  angles  resulting  from  quantum 
mechanical  treatment  of  scattering  in  the  lattice  potentials  in  the  secondary 
emission  medium. 

True  secondaries  are  electrons  that  are  emitted  from  the  conduction  or  va¬ 
lence  bands  of  the  atoms  comprising  the  impact  surface.  The  emitted  electron 
population  contains  about  90%  true  secondaries.  Energy  is  imparted  to  the 
lattice  electrons  over  some  time  long  compared  to  the  elastic  collision  time,  so 
the  electron  energy  distribution  can  be  modeled  as  a  Maxwell-Boltzmann  dis¬ 
tribution.  The  electrons  are  emitted  with  energies  from  a  thermal  distribution 
of  temperature  T: 


where  ks  is  Boltzmann’s  constant.  Due  to  the  timescale  of  the  emission  process, 
the  angle  of  emission  can  be  taken  to  be  isotropic: 

(5) 


1.3  Discussion 

The  secondary  model  described  above  has  been  implemented  in  part  in  the 
XPDP1  code  [7],  and  in  full  in  the  XOOPIC  code  [1].  In  the  scheme  described 
above,  the  secondary  emission  coefficient  can  indicate  that  a  fractional  number 
of  electrons  must  be  emitted.  The  above  codes  emit  fractional  electron  yields 
statistically,  using  a  random  number.  Another  technique  is  to  accumulate  frac¬ 
tional  particles  until  a  sufficient  level  is  achieved  to  emit  a  particle.  Variable 
particle  weighting  can  also  be  used  to  represent  fractional  electron  emission  [8]. 
The  secondary  electrons  are  typically  computed  after  each  electron  is  absorbed 
at  a  surface,  at  the  end  of  a  timestep,  and  advanced  into  the  simulation  in 
the  next  timestep  (XPDP1).  Riyopoulus  [9]  has  noted  a  banding  effect  when 
emitting  at  the  end  of  a  timestep,  due  to  an  average  emission  delay  of  At/2. 
This  was  repaired  by  computing  a  fractional  step  for  emission  which  accounts 
for  the  time  within  the  timestep  that  the  particle  impacted  the  emission  sur¬ 
face.  XOOPIC  is  capable  of  advancing  each  particle  the  remaining  fraction  of 
the  timestep,  eliminating  this  problem.  The  banding  effect  is  observed  only 
for  large  timesteps,  u)rfAt/2  &  1,  where  urf  is  the  driving  frequency  of  the 
multipactor  in  the  Riyopoulus  study. 
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Abstract 

A  digital  filtering  algorithm  for  particle  simulation  in  cylindrical  coordinates  is  derived. 
For  application  of  the  filter  to  the  charge  density,  the  filter  has  the  desirable  property  that 
total  charge  is  conserved,  and  the  filter  preserves  a  uniform  density.  Analysis  indicates  the 
filter  strength  decreases  with  radius,  and  that  hundreds  of  applications  are  required  to  reduce 
statistical  fluctuations  caused  by  finite  size  particles  weighted  to  the  mesh.  A  scheme  is  presented 
to  apply  n  passes  of  the  filter  using  a  single  matrix  multiplication  and  a  precomputed  matrix. 
The  filter  is  demonstrated  on  a  high  density  lamp  discharge  plasma  in  cylindrical  coordinates. 


1  Introduction 

Statistical  noise  has  long  been  a  significant  limitation  for  particle-in-cell  (PIG)  calculations  at  high 
densities.  This  noise  is  manifested  in  the  charge  and  current  density  source  terms  to  the  field 
equations,  which  results  in  fluctuations  in  fields.  Such  fluctuations  in  fields  can  cause  non-physical 
net  heating  of  particles  [1].  The  noise  and  consequent  fluctuations  are  inversely  proportional  to  the 
square  root  of  the  number  of  particles,  so  the  brute  force  reduction  of  noise  by  increasing  the  number 
of  particles  is  computationally  inefficient.  The  problem  is  significantly  more  severe  in  curvilinear 
coordinates,  where  the  particle  density  is  proportional  to  the  radius  r  in  cylindrical  coordinates, 
and  proportional  to  r 2  in  spherical  coordinates. 

A  common  remedy  for  this  problem  is  the  use  of  spectral  filters  or  digital  filters  [1],[2].  Spectral 
filters  operate  in  k-space,  most  appropriate  when  a  spectral  field  solution  is  performed.  Digital 
filters  are  more  amenable  to  finite  difference  field  solutions  in  real  space,  which  is  generally  the  case 
for  bounded  PIC  codes  in  curvilinear  coordinates. 

The  objective  of  this  work  is  to  extend  digital  filtering  schemes  to  cylindrical  coordinates.  It 
is  sufficient  to  consider  the  one-dimensional  cylindrical  case  here,  as  the  work  can  be  extended 
to  multiple  dimensions  in  a  straightforward  fashion  with  the  usual  algebraic  complications.  In 
Section  2,  the  algorithm  is  derived  and  an  implementation  is  described.  In  Section  3,  a  scheme  for 
applying  the  equivalent  of  multiple  passes  of  the  filter  with  a  single  matrix-vector  multiply  using 
a  pre-computed  filter  matrix  is  described.  In  Section  4,  the  algorithm  is  analyzed  for  both  single 
and  multiple  passes,  and  the  effect  of  the  filter  is  demonstrated.  In  Section  5,  the  conclusions  are 
presented. 
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2  Algorithm 

We  would  like  to  consider  a  digital  filter  in  cylindrical  coordinates  similar  to  the  1-2-1  spatial  filter 
in  Cartesian  coordinates  discussed  in  Birdsall  and  Langdon  [l] .  This  type  of  filter  is  useful  for 
reducing  the  short  wavelength  fluctuations  generated  by  discrete  particle  shapes  which  can  lead  to 
particle  heating.  Note  that  the  particle  noise  scales  as  the  square  root  of  the  number  of  particles; 
a  short  wavelength  filter  can  reduce  the  number  of  particles  required  to  meet  some  noise-induced 
heating  criterion  in  a  particle  simulation,  relaxing  demands  on  computing  time  to  achieve  the 
desired  noise  level. 

Consider  a  general  one-dimensional  mesh,  in  which  the  volumes,  V),  are  specified  for  0  <  j  <  N. 
For  physical  reasons,  we  require  conservation  of  charge,  Q =  Qn,  where  the  superscript  refers 
to  the  iteration  of  smoothing,  and 

=  (i) 

J-O 

We  must  choose  some  other  criterion  for  the  coefficients  of  the  filter.  For  a  cylindrical  system, 
it  is  possible  that  a  Bessel  function  might  be  a  more  natural  state,  but  here  for  simplicity  we  choose 
a  function  which  will  not  perturb  a  uniform  density  in  analogy  to  the  planar  case: 

Po+1  =  (!  -  «o)  Po  +  ctQPi, 

P]+1  =  aJPl- 1  +  (1  “  2 aj)  p"  +  otjpP+1,  0  <  j  <  N,  (2) 

P^+ 1  =  aNpnN_ l  +  (1  -  aN)  p%. 

Here,  aj  are  coefficients  to  be  determined  by  satisfying  charge  conservation.  We  can  see  by  inspec¬ 
tion  that  a  spatially  uniform  density  will  remain  uniform  after  filtering  by  Eq.  2. 

Applying  charge  conservation  to  Eq.  2,  we  obtain 

N  N 

E  viPnj  =  Y,ViPTl  =  M^-^)pno+aQp\]  +  VN[aNpnN_^{l-aN)pl} 
j— o  j"  o 

N-l 

+  T.vi 

j= 1 

Comparing  like  terms  of  we  obtain  a  set  of  linear  equations  for  the  coefficients: 


[Pj-i  +  Phi)  +  (!  “  2ai)  Pi 


(3) 


—Vqcxo  +  Viai  =  0, 

Vj-ictj-i  —  2  VjOLj  +  —  0, 

Vn-iqln-i  -  VNajy  =  0. 


(4) 


However,  these  are  not  linearly  independent,  as  consideration  of  the  3-row  case  demonstrates: 
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Clearly  the  middle  row  is  the  opposite  of  the  sum  of  the  other  two  rows.  This  means  we  must 
obtain  one  additional  equation.  If  we  place  a  constraint  on  the  filter  that  it  should  not  change 
the  sign  of  the  contribution  of  a  given  charge  density,  p™  to  the  co-located  filtered  density,  p^+1, 
then  we  constrain  0  <  oq  <  1.  So  that  this  algorithm  returns  the  familiar  (and  well-studied)  1-2-1 
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Figure  1:  Multipass  digital  filter  coefficients  for  a  100  pass  filter  on  a  20  by  20  uniformly  spaced 
radial  mesh  in  cylindrical  coordinates,  with  ao  =  1/2.  The  filter  coefficients  for  a  fixed  row 
(horizontal  slice)  show  the  contribution  of  the  unfiltered  data  to  a  particular  filtered  element, 
while  the  coefficients  for  a  fixed  column  (vertical  slice)  show  the  redistribution  profile  for  a  unit 
charge  at  the  appropriate  radius. 

digital  filter  for  uniform  meshes  in  Cartesian  coordinates,  we  might  choose  ao  =  1/2.  This  results 
in  the  generating  sequence  for  the  coefficients: 

ao  =  1/2, 

«i  =  aoVo/Vi, 

otj  =  2aj-iVj-i/Vj  -  aj-M-i/Vj,  2  <  j  <  N,  {  1 

CUN  =  ®N-  iVn-i/Vn- 

The  recursion  given  in  Eq.  6  yields  the  familiar  1-2-1  digital  filter  when  the  volumes  are  equal 
(except  the  end  cells,  which  have  half  the  volume  of  the  center  cells).  It  conserves  charge  for 
arbitrary  volumes,  and  preserves  a  uniform  distribution  for  arbitrary  volumes.  This  algorithm 
should  work  for  both  uniform  and  non-uniform  meshes  in  any  coordinate  system  in  one  dimension. 

3  Multiple  Applications  of  the  Digital  Filter 

As  shown  above,  the  particle  shape  is  weakly  impacted  by  the  smoothing  in  cylindrical  coordinates. 
Empirical  evidence  also  indicates  that  hundreds  of  passes  of  smoothing  are  required  to  reduce  noise 
in  discharges  in  order  to  reduce  numerical  heating.  Since  the  algorithm  can  become  expensive  with 
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hundreds  of  passes  of  Eq.  2,  we  consider  a  scheme  by  which  n  passes  of  the  filter  can  be  applied 
simultaneously.  Taking  the  filter  to  be  of  the  form 

pn+1  =  Ap\  (7) 


where  pn  is  the  charge  density  vector  over  the  spatial  dimension  after  n  applications  of  the  filter, 
and  A  is  the  matrix  constructed  from  aj  given  by 


A  = 
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&N- 1 
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(8) 


When  we  apply  multiple  passes  of  the  digital  filter,  we  simply  nest  applications  of  the  linear 
operation  given  in  Eq.  7.  This  leads  to  the  easily  implemented  and  efficient  multiple  pass  algorithm, 

Pn  =  Anp° ,  (9) 


where  An  can  be  precomputed  for  applications  where  the  digital  filter  level  is  prescribed  initially. 

The  coefficient  matrix  for  a  20  by  20  uniformly  spaced  mesh,  with  the  initial  coefficient  an  =  1/2, 
is  shown  in  Fig.  1.  Slices  with  fixed  row  numbers  characterize  the  radial  spread  of  the  filtered 
quantity  at  a  given  radius,  from  which  we  see  that  the  filtering  produces  broadening  of  an  initial 
distribution  which  increases  with  decreasing  radius.  Slices  with  fixed  column  numbers  show  the 
contributions  of  the  pre-filtered  components  to  a  component  at  the  radius  corresponding  to  the 
column. 


4  Discussion 

In  planar  coordinates,  the  shape  of  a  linearly  weighted  charge  in  one  dimension  remains  linear, 
but  the  width  of  the  charge  increases  in  proportion  to  the  decrease  in  the  peak.  In  cylindrical 
coordinates,  the  analogous  weighting,  found  from  recursing  Eq.  6  with  a:o  =  1/2,  results  in  a  much 
milder  impact  on  the  particle  shape.  The  particle  shapes  for  the  initial  linear  spline  and  1,  10, 
and  100  applications  of  the  an  =  1/2  filter  are  shown  in  Fig.  2.  A  single  application  of  the  filter 
is  almost  indistinguishable  from  the  original  linear  particle  shape.  After  10  applications  a  small 
change  in  shape  is  discernible,  and  after  100  applications  the  shape  is  considerably  smoother.  Note 
that  the  filtered  shapes  have  a  discontinuity  in  the  first  derivative  at  each  cell  edge,  and  the  second 
derivative  on  the  left  side  of  the  particle  is  slightly  positive,  while  it  is  slightly  negative  on  the  right 
side. 

Spectral  analysis  corresponding  to  the  above  figure. 

The  multi-pass  digital  filter  has  been  applied  to  a  high  density  1  cm  cylindrical  discharge  using 
the  code  XPDCl.  The  discharge  is  a  model  of  a  fluorescent  lamp  discharge,  which  suffers  from  sig¬ 
nificant  numerical  heating  due  to  the  statistical  fluctuations  which  are  approximately  proportional 
to  1/r.  This  unfavorable  radial  scaling  makes  the  heating  problem  very  expensive  to  address  by 
increasing  the  number  of  particles,  Np,  since  the  noise  level  is  proportional  to  .  The  charge 

density  from  the  multi-pass  digital  filter  with  n  =  400  is  compared  to  the  initial  charge  density 
obtained  from  linear  weighting  in  Fig.  3.  Note  the  1/r  characteristic  envelope  on  the  unfiltered 
charge  density,  and  the  reduction  in  short  wavelength  fluctuations  with  the  filter. 
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Figure  2:  Shape  of  a  standard  linearly  weighted  particle,  and  shape  after  1,  10,  and  100  passes  of 
the  qiq  =  1/2  cylindrical  filter.  The  linear  spline  and  1-pass  curve  are  nearly  co-located. 

5  Conclusions 

A  digital  filter  for  discrete  functions  such  as  charge  density  was  derived  for  cylindrical  coordinates. 
The  filter  is  designed  to  have  the  desirable  characteristics  of  charge  conservation,  as  well  as  pre¬ 
serving  a  uniform  density.  Although  single  passes  of  the  filter  have  only  a  small  effect,  a  scheme 
is  described  which  can  apply  the  equivalent  of  an  arbitrary  number  of  filter  passes  with  a  single 
matrix- vector  multiplication  using  a  sparse  matrix  computed  initially.  The  digital  filter  has  a  small 
per-pass  effect,  but  can  have  a  significant  multi-pass  effect.  The  digital  filter  is  demonstrated  to  be 
effective  in  reducing  statistical  fluctuations  in  the  charge  density  for  a  high  density  lamp  discharge 
which  would  suffer  from  numerical  heating  without  the  filter. 
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Charge  Density (r)  (C/mA-3) 


Figure  3:  Charge  density  from  an  n  =  400  multipass  filter  compared  to  unfiltered  data.  The 
discharge  is  net  neutral  except  in  the  sheath,  so  all  short  wavelength  fluctuations  inside  the  sheath 
radius  are  attributed  to  statistical  noise. 
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Abstract 

A  general  method  for  computing  charge  and  current  density  source  terms  for  Maxwell’s  equa¬ 
tions  from  particles  weighted  to  a  mesh  is  described.  The  method  presented  here  eliminates  the 
need  for  correction  factors  often  applied  in  curvilinear  coordinates  to  compensate  for  errors  at 
the  edge  of  the  system,  and  in  the  interior  as  well  for  nonuniform  meshes.  Generality  is  achieved 
by  weighting  volume  elements  using  a  spline  symmetric  to  that  by  which  particle  charge  and 
current  are  weighted  to  the  mesh.  The  method  presented  has  a  number  of  desirable  proper¬ 
ties,  including  conservation  of  charge,  preservation  of  a  uniform  distribution,  and  generality  on 
nonuniform  meshes  with  arbitrary  particle-mesh  interpolation  schemes.  The  method  recovers 
the  exact  answer  in  the  limit  of  mesh  sizes  approaching  zero. 

Keywords:  PIC,  particle,  weighting,  radial,  cylindrical,  spherical,  correction 

1  Introduction 

In  charge  and  current  accumulation  schemes  commonly  used  in  particle-in-cell  codes,  a  systematic 
error  occurs  on  the  boundary  cells  in  curvilinear  coordinates.  The  error  is  particularly  severe  on 
the  axis  in  cylindrical  and  spherical  models,  with  systematic  errors  of  33%  and  100%  larger  than 
the  theoretical  value  for  either  charge  or  current  density  in  cylindrical  and  spherical  coordinates, 
respectively  Additional  error  occurs  at  the  outer  edge  of  the  system,  and  throughout  the  interior 
for  non-uniform  meshes.  This  error  leads  to  similar  errors  in  the  forces  calculated  at  those  points, 
as  well  as  exacerbating  radial  noise  due  to  particle  statistics,  which  is  largest  near  r  =  0.  The  errors 
in  density  do  not  decrease  with  decreasing  mesh  size.  This  problem  is  not  present  in  Cartesian 
coordinates,  even  for  non-uniform  meshes. 
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The  issue  of  radial  correction  factors  for  the  charge  density  term  has  been  addressed  in  recent 
works  by  Ruyten  [l]  and  by  Larsen  et  Ah.  [2].  These  works  derive  the  error  for  linear  weighting  to 
uniform  meshes  in  cylindrical  coordinates  and  propose  correction  terms.  In  [1],  the  author  states 
that  the  axial  behavior  is  beyond  the  scope  of  the  work.  In  [2],  the  analysis  is  carried  out  including 
the  origin,  and  the  proposed  correction  is  shown  to  approach  the  theoretical  value  for  large  numbers 
of  particles.  We  will  demonstrate  that  the  result  of  [2]  is  exact  for  the  uniform  charge  density  case 
despite  errors  shown  in  that  paper  for  charge  density  at  the  axis. 

In  this  work,  the  error  is  analyzed  for  the  uniform  particle  density  distribution  for  the  charge 
density,  p.  A  general  scheme  is  developed  which  applies  the  same  interpolation  scheme  to  both 
the  particle  and  the  volume  or  surface  element  in  order  to  obtain  a  density  which  eliminates  the 
systematic  error  in  curvilinear  coordinates.  The  new  scheme  is  compared  with  analytic  theory  in 
cylindrical  and  spherical  coordinates.  The  algorithm  is  extended  to  the  current  density,  J. 


2  Density  Errors 


In  this  section,  the  errors  in  the  most  common  interpolation  scheme  used  in  particle  simulation  [3], 
linear  weighting,  are  demonstrated  for  uniform  density  in  one  dimensional  on  a  nonuniform  mesh 
in  cylindrical  and  spherical  coordinates.  For  an  arbitrary  continuum  particle  distribution  specified 
by  /(?’),  the  density  in  cylindrical  coordinates  is  given  by: 

f(r')drf 


n(r)  —  — 


•r~\-dr 
r 


(i) 


jr+dr  2irrt(lri 

For  the  uniform  particle  distribution,  f(r)  =  27rr,  and  we  obtain  n(r )  =  1.  Similarly,  in  spherical 
coordinates  for  the  uniform  particle  distribution  f(r)  —  47rr2  we  obtain: 


Jrr+dr  f(r')dr’ 


n  (r)  = 


=  1 
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jr+dr  ^ 

For  linear  weighting  on  a  nonuniform  mesh  in  cylindrical  coordinates,  we  can  write  the  density 
at  an  arbitrary  intermediate  node  for  an  arbitrary  particle  distribution  f(r)  in  the  standard  way 
[3]: 

_  1  f(r)y-^-^dr  +  J? hl  /(r)gf 


frj+ 1/2  27t rdr 


(3) 


where  the  rj  refer  to  the  position  of  the  jth  mesh,  and  rJ+i/2  =  (r,-  +  r^+i)  /2.  The  edge  densities 
are  written  (valid  also  on  axis,  ro  =  0): 
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For  a  uniform  particle  distribution,  f(r)  =  27rr,  the  standard  discrete  densities  in  Eqs.  3-5 
become: 

4  r3~  1  +  r3  +  TH1 
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On  a  non-uniform  mesh,  the  linear  weighting  with  uncorrected  volumes  always  produces  the  incor¬ 
rect  result  for  all  cells.  For  a  uniform  mesh,  we  see  from  inspection  of  Eq.  6  that  the  correct  solution 
is  produced  for  the  interior  cells,  but  the  systematic  error  persists  at  the  edges.  For  nonuniform 
meshes,  the  error  occurs  at  interior  points  as  well.  The  systematic  error  is  1/3  on  axis  for  meshes 
which  include  the  axis,  and  is  independent  of  grid  spacing  Ar.  Note  that  the  error  in  the  outer 
edge  is  small  for  N  1.  The  uncorrected  results  for  a  uniform  and  nonuniform  mesh  are  plotted 
in  Fig.  1. 

A  similar  development  for  linear  weighting  in  spherical  coordinates  for  the  uniform  particle 
distribution,  f(r)  =  47rr2,  gives: 


n * 


=  2  r)_ i  +  rj-xr,-  +  r^rj+i  +  r)  +  rjrj+1  +  r?+1 
rj_l  +  Zrj-irj  +  rj-irj+i  +  3 r)  +  3 •rjrJ+1  +  rj+1 : 

,3  rl  +  2r0n  +  r\ 


no  =  2 
n-N  =  2  2 


„  ■>  ,  - o,  and 

7tq  +  4r0ri  +  r\ 

rAr-l  +  2rjy_ir  at  +  3  r% 
_l  +  4rAr_arAr  +  7r%’ 


N 


For  uniform  mesh  spacing,  Ar,  Eqs.  9-11  become: 

6  r?  +  Al 


Uj  2 12r?  +  A?  ’ 


n0 


2ftj±»oAL±A^,  ond 


n^v  —  2 


12r^  +  6r0Ar  +  A^  ’ 
6r^  -  4r^Ar  +  A2r 


12ri  -  6rNAr  +  A% ' 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 


Note  that  the  uncorrected  density  for  the  spherical  case  is  always  in  error,  approaching  the  correct 
solution  only  for  rj  »  Ar.  The  uncorrected  results  for  a  uniform  and  nonuniform  mesh  are  plotted 
in  Fig.  2.  The  systematic  error  is  100%  on  axis  for  meshes  which  include  the  axis,  independent  of 
grid  spacing  Ar. 
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3  Improved  Density  Weighting  Algorithm 


In  this  section,  a  general  algorithm  is  developed  for  computing  the  charge  density  in  non-uniform 
meshes  for  curvilinear  coordinates  using  arbitrary  particle-mesh  interpolation.  Consider  a  modified 
calculation  for  the  volume  in  which  differential  volume  elements  are  weighted  to  the  mesh  using 
the  same  algorithm  as  the  charge  weighting.  A  general  method  is  suggested  by  the  definition  of 
density  in  cylindrical  (Eq.  1)  and  spherical  (Eq.  2)  coordinates,  weighting  both  the  charge  and  the 
volume: 

.  Jr  fJ^Wjjyjdr 

j  St  Wj  (r)  dV  ’  1  j 

where  Wj  (r)  is  an  interpolation  function  which  weights  particles  at  position  r  to  mesh  j,  and  dV  is 
a  volume  element,  given  in  one  dimension  by  dV  =  2ttt  in  cylindrical  coordinates  and  dV  —  47 it2 
in  spherical  coordinates. 

For  linear  weighting  in  cylindrical  coordinates,  the  density  can  then  be  written: 

■J?-,  *  M%g;dr 


Uj  = 


f^r^dr  +  f^lvr^dr 
The  edge  densities  are  found  simply  by  dropping  the  out  of  bounds  integrals: 


(16) 


no  = 


rri 

Jrp  J  V  /  ri  —  rp 

O  ‘I'KT^dr  ’ 

Jrn  r i  ~~rn 


and 


UN 


frN  f(r)  — 
JrN-l  J  \  /  Tat  — 


r-rN- 1 


rN-rN-i 


dr 


W 2nr 


dr 


(17) 


(18) 


VN-TN-1 

For  /(r)  =  27rr,  we  obtain  the  exact  solution  for  all  0  <  j  <  N,  rtj  =  1.  For  more  general 
particle  distributions  (specifically  when  /  is  not  a  linear  function  of  r),  Eqs.  16-18  result  in  an 
error  term  which  depends  upon  the  weighting  scheme;  for  example  the  linear  weighting  function 
results  in  a  error  proportional  to  A^. 

This  can  be  implemented  in  a  standard  linear  weighting  scheme  simply  by  modifying  the  volumes 
used  to  compute  rij.  Taking  the  cell  size  in  the  2-direction  to  be  Az,  the  volumes  become: 


7 r 


Vj  =  Az—  [tj+i  (rj  +  rj+1)  -  rj- 1  (r,-_ i  +  rj)} , 


(19) 


7T 


V0  =  Az-  (rj  -  r0)  (2 r0  +  n) ,  and 


7T 


Vn  =  Az—  (rjv  —  vn-i)  (r^-i  +  2 r^) . 
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(20) 

(21) 
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The  corrected  and  uncorrected  densities  have  been  implemented  in  the  XOOPIG  [4]  and  XPDC1 
[5]  codes,  and  comparisons  of  the  results  are  shown  in  Fig.  1.  The  comparisons  shown  are  valid 
when  there  is  at  least  1  particle  per  cell.  Total  charge  and  volume  are  also  identically  conserved 
with  the  corrected  method.  This  method  can  be  easily  extended  to  arbitrary  weighting  functions 
by  using  the  desired  weighting  in  Eqs.  16-18. 

For  linear  weighting  in  spherical  coordinates  the  density  can  be  written: 


IrLi  4?rr2£Sdr  +  Ir-+1  ^7^-— dr 


rj+1-rj 


with  edge  densities  again  obtained  by  dropping  out  of  bounds  integrals. 

For  the  uniform  particle  distribution,  /(r)  =  47rr2,  we  again  obtain  the  exact  solution  for  all 
0  <  j  <  AT,  rij  —  1.  For  more  general  particle  distributions  we  obtain  the  same  result  as  that  of  the 
cylindrical  case  above. 


Similar  to  the  cylindrical  case,  the  method  can  be  implemented  for  the  spherical  scheme  by 
precomputing  the  volumes  using  the  denominator  of  Eq.  22: 


Vj  =  J  (ri+i  -  rj- 1)  (rj- 1  +  rj-irj  +  rj-iTj+1  +  rf  +  rjrj+1  +  rj+1)  ,  (23) 

V0  =  ^  (ri  -  r0)  (Sr'o  +  2rxr0  +  rfj  ,  and  (24) 

Vn  —  g-  (xn  —  rN-i)  (rff-i  +  ^n-i^n  +  3 .  (25) 

The  densities  computed  with  the  method  presented  here  are  compared  to  the  uncorrected  den¬ 
sities  for  linear  weighting  on  a  uniform  and  nonuniform  mesh  in  Fig.  2. 


4  Current  Density 

This  algorithm  can  be  applied  to  current  density  for  the  electromagnetic  source  term  in  a  straight¬ 
forward  manner  by  using  J  =  nqv.  The  quantity  qv  is  weighted  to  the  mesh  for  each  particle,  and 
the  result  is  divided  by  the  surface  area,  S.  To  obtain  the  surface  area,  consider  the  differential 
elements  of  the  components  in  cylindrical  coordinates, 

dSr  =  rdddz ,  (26) 


dSg  =  clrdz,  and 

(27) 

dSz  =  rdrdQ. 

(28) 
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It  is  evident  from  Eqs.  26-28  that  only  the  axial  component  of  S  will  result  in  a  nonlinear 
dependence  on  r  when  computing  the  surface  area  from  S  =  fW dS ,  where  W  is  the  particle-mesh 
weighting  used.  The  results  for  the  axial  component  of  the  current  density  are  trivially  different 
from  the  charge  density  results  derived  above,  and  the  comparison  plots  are  identical. 

5  Conclusions 

An  algorithm  for  obtaining  the  correct  charge  and  current  densities  in  curvilinear  coordinate  systems 
for  arbitrary  particle  interpolation  schemes  is  described.  Volumes  and  surface  areas  are  weighted  to 
the  mesh  using  the  same  interpolation  scheme  used  to  weight  particles.  The  method  recovers  the 
charge  density  correction  factors  for  linear  weighting  in  cylindrical  coordinates  [2] ,  but  extends  to 
the  more  general  case.  The  method  has  the  notable  beneficial  properties  of  conserving  charge  and 
current  as  well  as  total  volume  and  surface  area  on  a  general  orthogonal  mesh.  The  algorithm  yields 
the  exact  answer  for  the  uniform  particle  distribution,  as  demonstrated  here;  arbitrary  distributions 
are  similarly  correct  to  the  mesh  resolution,  becoming  exact  as  the  mesh  size  appi'oaches  zero  for 
continuum  distributions.  This  is  a  significant  departure  of  the  pervious  scheme  of  computing  using 
incorrect  volumes  and  then  using  correction  factors  designed  for  special  cases  [2].  The  results  are 
easily  extended  to  two  and  three  dimensions. 
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Figure  1:  Mesh  densities  computed  using  standard  linear  weighting  (uncorrected)  on  a  uni¬ 
form  and  nonuniform  mesh  in  cylindrical  coordinates,  as  well  as  the  corrected  density  on  both 
meshes.  The  uniform  mesh  has  TV  =  10  steps;  the  nonuniform  mesh  is  computed  from  Vj  = 
+  (rN  —  ro)  such  that  the  interval  contains  TV  =  10  equal  volumes. 
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Figure  2:  Mesh  densities  computed  using  standard  linear  weighting  (uncorrected)  on  a  uniform  and 

nonuniform  mesh  in  spherical  coordinates,  as  well  as  the  corrected  density  on  both  meshes.  The  uni- 

r  -i  1/3 

form  mesh  has  N  =  10  steps;  the  nonuniform  mesh  is  computed  from  Tj  =  +  (r^  —Tq)/N  , 

such  that  the  interval  contains  N  =  10  equal  volumes. 
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A  general  method  for  computing  charge  and  current  density  source  terms 
for  Maxwell’s  equations  from  particles  weighted  to  a  mesh  is  described.  The 
method  presented  here  eliminates  the  need  for  correction  factors  often  ap¬ 
plied  in  curvilinear  coordinates  to  compensate  for  errors  at  the  edge  of 
the  system,  and  in  the  interior  as  well  for  nonuniform  meshes.  Gener¬ 
ality  is  achieved  by  weighting  volume  elements  using  a  spline  symmetric 
to  that  by  which  particle  charge  and  current  are  weighted  to  the  mesh. 
The  method  presented  has  a  number  of  desirable  properties,  including  con¬ 
servation  of  charge,  preservation  of  a  uniform  distribution,  and  generality 
on  nonuniform  meshes  with  arbitrary  particle- mesh  interpolation  schemes. 
The  method  recovers  the  exact  answer  in  the  limit  of  mesh  sizes  approach¬ 
ing  zero. 
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1.  INTRODUCTION 

In  charge  and  current  accumulation  schemes  commonly  used  in  particle-in-cell 
codes,  a  systematic  error  occurs  on  the  boundary  cells  in  curvilinear  coordinates. 
The  error  is  particularly  severe  on  the  axis  in  cylindrical  and  spherical  models,  with 
systematic  errors  of  33%  and  100%  larger  than  the  theoretical  value  for  either  charge 
or  current  density  in  cylindrical  and  spherical  coordinates,  respectively.  Additional 
error  occurs  at  the  outer  edge  of  the  system,  and  throughout  the  interior  for  non- 
uniform  meshes.  This  error  leads  to  similar  errors  in  the  forces  calculated  at  those 
points,  as  well  as  exacerbating  radial  noise  due  to  particle  statistics,  which  is  largest 
near  r  =  0.  The  errors  in  density  do  not  decrease  with  decreasing  mesh  size.  This 
problem  is  not  present  in  Cartesian  coordinates,  even  for  non-uniform  meshes. 

The  issue  of  radial  correction  factors  for  the  charge  density  term  has  been  ad¬ 
dressed  in  recent  works  by  Ruyten  [1]  and  by  Larsen  et  al.  [2].  These  works  derive 
the  error  for  linear  weighting  to  uniform  meshes  in  cylindrical  coordinates  and  pro¬ 
pose  correction  terms.  In  [1],  the  author  states  that  the  axial  behavior  is  beyond 
the  scope  of  the  work.  In  [2],  the  analysis  is  carried  out  including  the  origin,  and  the 
proposed  correction  is  shown  to  approach  the  theoretical  value  for  large  numbers  of 
particles.  We  will  demonstrate  that  the  result  of  [2]  is  exact  for  the  uniform  charge 
density  case  despite  errors  shown  in  that  paper  for  charge  density  at  the  axis. 

In  this  work,  the  error  is  analyzed  for  the  uniform  particle  density  distribution 
for  the  charge  density,  p.  A  general  scheme  is  developed  which  applies  the  same 
interpolation  scheme  to  both  the  particle  and  the  volume  or  surface  element  in 
order  to  obtain  a  density  which  eliminates  the  systematic  error  in  curvilinear  co¬ 
ordinates.  The  new  scheme  is  compared  with  analytic  theory  in  cylindrical  and 
spherical  coordinates.  The  algorithm  is  extended  to  the  current  density,  J. 


2.  DENSITY  ERRORS 
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In  this  section,  the  errors  in  the  most  common  interpolation  scheme  used  in 
particle  simulation  [3],  linear  weighting,  are  demonstrated  for  uniform  density  in 
one  dimensional  on  a  nonuniform  mesh  in  cylindrical  and  spherical  coordinates.  For 
an  arbitrary  continuum  particle  distribution  specified  by  /(r),  the  exact  density  in 
cylindrical  coordinates  is  given  by: 


n(r) 


/;+dr/(rQrfr' 
Jrr+rfr  27 rr'dr' 


(1) 


For  the  uniform  particle  distribution,  f(r)  =  2ir r,  and  we  obtain  n(r)  =  1. 
Similarly,  in  spherical  coordinates  for  the  uniform  particle  distribution  f(r)  =  47rr2 
we  obtain: 


n  (r)  = 


£+drmdr' 
fr+dr  47t  (r')2  dr1 


-  1 


(2) 


In  the  classical  particle  scheme  [3] ,  the  cell  volumes  are  computed  from  geometric 
considerations  independently  of  the  weighting  scheme  used  to  accumulate  charge 
to  the  grid.  Because  the  volume  elements  are  nonlinear  in  the  spatial  variable  for 
curvilinear  coordinates,  this  leads  to  systematic  errors.  This  conceptual  error  is 
precisely  why  the  previous  works  ([1]  and  [2])  describe  corrections  to  the  weighting 
scheme.  For  linear  weighting  on  a  nonuniform  mesh  in  cylindrical  coordinates,  we 
can  write  the  density  at  an  arbitrary  intermediate  node  for  an  arbitrary  particle 
distribution  f(r)  in  the  standard  way  [3]  in  order  to  make  the  conceptual  error 
more  apparent: 


frj  f(r)  r  r*  1  dr  +  frj+1  f(r)  r**1  -^-dr 

—  ZlLzl  J  V  ^  r j  ~  ri  - 1  JrJ  _  V  /ri+i~rj 
2yrr<^r 


(3) 


where  the  Vj  refer  to  the  position  of  the  jth.  mesh,  and  T/+1/2  =  (rj  4*  /2.  The 

edge  densities  are  written  (valid  also  on  axis,  Tq  =  0): 


no  = 


fri  /( r — —  dr 
Jr 0  J  \  /  ri-rp 

fri/2  27t rdr 

Jr  0 


and 


(4) 
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nN  = 


[  /(r)-2^ 

Jrjv-i  J  v  1  r  at  —  r 


-dr 


frN  2irrdr 

JT  pf  —  1/2 


(5) 


For  a  uniform  particle  distribution,  /(r)  =  27rr,  the  standard  discrete  densities 


in  Eqs.  3-5  become: 


4  rj-i  +  rj  +rj+i 

3  7'j  —  i  +  2  Tj  +Tj+ 1  ’ 

(6) 

4  2r0  +  rx 

n0  —  - - ,  and 

3  3ro  +  ri 

(7) 

4rN-i  +  2rN 
nN~  3rw_1  +  3rjV- 

(8) 

On  a  non-uniform  mesh,  the  linear  weighting  with  uncorrected  volumes  always  pro¬ 
duces  the  incorrect  result  for  all  cells.  For  a  uniform  mesh,  we  see  from  inspection 
of  Eq.  6  that  the  correct  solution  is  produced  for  the  interior  cells,  but  the  system¬ 
atic  error  persists  at  the  edges.  For  nonuniform  meshes,  the  error  occurs  at  interior 
points  as  well.  The  systematic  error  is  1/3  on  axis  for  meshes  which  include  the 
axis,  and  is  independent  of  grid  spacing  Ar.  Note  that  the  error  in  the  outer  edge 
is  small  for  TV  »  1.  The  uncorrected  results  for  a  uniform  and  nonuniform  mesh 
are  plotted  in  Fig.  1. 

A  similar  development  for  linear  weighting  in  spherical  coordinates  for  the  uni¬ 
form  particle  distribution,  f(r)  =  4-zrr2,  gives: 

_  rj_  i  +  Tj-iTj  +  rj-iTj+ 1  +  rj  +  rjrj+i  +  rf+1 
r)_ ,  +  3  Tj-iTj  +  Tj-iTj+i  +  3  r?  +  3r;rJ+i  +  rj+1  ’ 


_  3rg  +  2r^rH-r| 
n°  7 r%  +  4r0ri  +  r\  ’ 


and 


_  0rw— l  +2rAr-i7>  +3rAr 
nN  ~  r%_1+4rN-irN +7r% 

For  uniform  mesh  spacing,  Ar,  Eqs.  9-11  become: 


„6r?  +  A2 

=  2T2^Ta 2 


(10) 

(11) 


(12) 


DRAFT 


August  13,  2001,  1:32pm  DRAFT 


TITLE  HEAD 


5 


_  6rg  +  4r0Ar  +  A? 
—  12  rl  +  6r0Ar  +  ’ 


and 


(13) 


_  0  Gr%  -  4rAr  A,.  +  A^ 

nN~Zl2r%-6rNAr  +  Ar 
Note  that  the  uncorrected  density  for  the  spherical  case  is  always  in  error,  approach¬ 
ing  the  correct  solution  only  for  Tj  >  Ar.  The  uncorrected  results  for  a  uniform 
and  nonuniform  mesh  are  plotted  in  Fig.  2.  The  systematic  error  is  100%  on  axis 
for  meshes  which  include  the  axis,  independent  of  grid  spacing  Ar. 


3.  IMPROVED  DENSITY  WEIGHTING  ALGORITHM 

In  this  section,  a  general  algorithm  is  developed  for  computing  the  charge  density 
in  non-uniform  meshes  for  curvilinear  coordinates  using  arbitrary  particle-mesh 
interpolation.  Consider  a  modified  calculation  for  the  volume  in  which  differential 
volume  elements  are  weighted  to  the  mesh  using  the  same  algorithm  as  the  charge 
weighting.  A  general  method  is  suggested  by  the  definition  of  density  in  cylindrical 
(Eq.  1)  and  spherical  (Eq.  2)  coordinates,  weighting  both  the  charge  and  the 
volume: 


n\ 


fr  f(r)Wj  (r)  dr 


(15) 


J  JrW}(r)dV  ’ 

where  Wj  (r)  is  an  interpolation  function  which  weights  particles  at  position  r  to 
mesh  j,  and  dV  is  a  volume  element,  given  in  one  dimension  by  dV  =  2irr  in 
cylindrical  coordinates  and  dV  =  47rr2  in  spherical  coordinates. 

For  linear  weighting  in  cylindrical  coordinates,  the  density  can  then  be  written: 


n7* 


n  fir)  r~ri-1  dr  +  fVj+1  f(r)d2±^Ldr 
iZlZllzl  V  r i + 1  ~ r J 

P  2tt rp^-dr  4-  P*+1  2irr^^-dr  ‘ 

Jrj-i  t  j  Tj  —  i  Jrj  rj+i-rj 


(16) 


The  edge  densities  are  found  simply  by  dropping  the  out  of  bounds  integrals: 


n0  = 


(17) 
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and 


JrN- /  I'N  —  rN-i  /1  o\ 

=  r  2„r  (18) 

Jr^- 1  rjv~riv-i 

For  f(r)  —  27rr,  we  obtain  the  exact  solution  for  all  0  <  j  <  TV,  rtj  —  1.  For 
more  general  particle  distributions  (specifically  when  /  is  not  a  linear  function  of 
7’) ,  Eqs.  16-18  result  in  an  error  term  which  depends  upon  the  weighting  scheme; 
for  example  the  linear  weighting  function  results  in  a  error  proportional  to  A2. 

This  can  be  implemented  in  a  standard  linear  weighting  scheme  simply  by  mod¬ 
ifying  the  volumes  used  to  compute  nj .  Taking  the  cell  size  in  the  ^-direction  to  be 
AZf  the  volumes  become: 


7 r 


Vj  =  Az—  [rj+1  {r3  +  rj+1)  -  rj-X  (ry  l  +rj)], 


(19) 


Vo  =  A*  -  (rx  -  r0)  (2r0  +  n) ,  and 


(20) 


V} v  =  Az—  (rpr  —  rjM-i)  (r^-i  +  2r^) . 


(21) 


The  corrected  and  uncorrected  densities  have  been  implemented  in  the  XOOPIC 
[4]  and  XPDC1  [5]  codes,  and  comparisons  of  the  results  are  shown  in  Fig,  1. 
The  comparisons  shown  are  valid  when  there  is  at  least  1  particle  per  cell.  Total 
charge  and  volume  are  also  identically  conserved  with  the  corrected  method.  This 
method  can  be  easily  extended  to  arbitrary  weighting  functions  by  using  the  desired 
weighting  in  Eqs.  16-18. 

For  linear  weighting  in  spherical  coordinates  the  density  can  be  written: 


ri4 


(22) 


3  f7j  4^7  *2JLJjLJ-dr  +  CJ+1  47rr2rrt-— r--rfr  ’ 

Jrj-i  Tj  Vj  —  i  Jrj  ri+i  rj 

with  edge  densities  again  obtained  by  dropping  out  of  bounds  integrals. 

For  the  uniform  particle  distribution,  /(r)  —  47rr2,  we  again  obtain  the  exact 
solution  for  all  0  <  j  <  TV,  n,j  —  1.  For  more  general  particle  distributions  we 
obtain  the  same  result  as  that  of  the  cylindrical  case  above. 
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Similar  to  the  cylindrical  case,  the  method  can  be  implemented  for  the  spherical 
scheme  by  precomputing  the  volumes  using  the  denominator  of  Eq.  22: 

Vj  =  |0j+1  -rv-i)  (rj_,  4-  rj-irj  +rj_,ri+1  +r]  +rjrj+1  +r|+1)  ,  (23) 

Vo  =  |  (r i  -  r0)  (3ro  +  2r!r0  4-  ?’j) ,  and  (24) 

vn-j  (rN  -  rN- 1)  (rN-i  +  ^N-\rN  +  3 r%)  .  (25) 

The  densities  computed  with  the  method  presented  here  are  compared  to  the 
uncorrected  densities  for  linear  weighting  on  a  uniform  and  nonuniform  mesh  in 
Fig.  2. 

4.  CURRENT  DENSITY 

This  algorithm  can  be  applied  to  current  density  for  the  electromagnetic  source 
term  in  a  straightforward  manner  by  using  J  =  nqv.  The  quantity  qv  is  weighted 
to  the  mesh  for  each  particle,  and  the  result  is  divided  by  the  surface  area,  S. 
To  obtain  the  surface  area,  consider  the  differential  elements  of  the  components  in 
cylindrical  coordinates, 


dSr  =  rdOdz , 

(26) 

dS$  =  drdz ,  and 

(27) 

dSz  —  rdrdO . 

(28) 

It  is  evident  from  Eqs.  26-28  that  only  the  axial  component  of  S  will  result  in 
a  nonlinear  dependence  on  r  when  computing  the  surface  area  from  S  =  f  WdS , 
where  W  is  the  particle-mesh  weighting  used.  The  results  for  the  axial  component 
of  the  current  density  are  trivially  different  from  the  charge  density  results  derived 
above,  and  the  comparison  plots  are  identical. 

5.  CONCLUSIONS 
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An  algorithm  for  obtaining  the  correct  charge  and  current  densities  in  curvilin¬ 
ear  coordinate  systems  for  arbitrary  particle  interpolation  schemes  is  described. 
Volumes  and  surface  areas  are  weighted  to  the  mesh  using  the  same  interpolation 
scheme  used  to  weight  particles.  The  method  recovers  the  charge  density  correc¬ 
tion  factors  for  linear  weighting  in  cylindrical  coordinates  [2],  but  extends  to  the 
more  general  case.  The  method  has  the  notable  beneficial  properties  of  conserving 
charge  and  current  as  well  as  total  volume  and  surface  area  on  a  general  orthogonal 
mesh.  The  algorithm  yields  the  exact  answer  for  the  uniform  particle  distribution, 
as  demonstrated  here.  Arbitrary  distributions  are  similarly  correct  to  the  mesh 
resolution,  becoming  exact  as  the  mesh  size  approaches  zero  for  continuum  distri¬ 
butions;  the  corresponding  electric  field  can  also  be  shown  to  approach  the  exact 
solution  in  the  same  limit.  This  is  a  significant  conceptual  departure  of  the  previ¬ 
ous  scheme  of  computing  using  incorrect  volumes  and  then  using  correction  factors 
designed  for  special  cases  [2].  The  results  are  easily  extended  to  two  and  three 
dimensions. 
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FIG.  1.  Mesh  densities  computed  using  standard  linear  weighting  (uncorrected)  on  a 


uniform  and  nonuniform  mesh  in  cylindrical  coordinates,  as  well  as  the  corrected  density  on 


both  meshes.  The  uniform  mesh  has  TV  =  10  steps;  the  nonuniform  mesh  is  computed  from 


rj  —  -f  (r^  —  rg)  /TV,  such  that  the  interval  contains  TV  =  10  equal  volumes. 
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FIG.  2.  Mesh  densities  computed  using  standard  linear  weighting  (uncorrected)  on  a 


uniform  and  nonuniform  mesh  in  spherical  coordinates,  as  well  as  the  corrected  density  on  both 
meshes.  The  uniform  mesh  has  jV  =  10  steps;  the  nonuniform  mesh  is  computed  from  rj  = 
(V|_  1  +  (r^-  —  r§)  /jvj  ^  ,  such  that  the  interval  contains  N  =  10  equal  volumes. 
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